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SUMMARY 


This  final  report  summarizes  the  work  performed  under  contract  No. 
N00014-84-C-0408  to  Office  of  Naval  Research  and  Is  accompanied  by  the 
software  listings  of  the  developed  code. 

The  goal  of  the  multipath  project  Is  to  develop  techniques  to  make  use  of 
the  multipath  relectlons  of  acoustic  signals  travelling  In  the  ocean  to 
localize  and  track  radiating  sources.  There  have  been  two  main  areas  of  study 
under  this  contract.  First  area  concentrates  on  the  detection  and  estimation 
of  the  multipath  (or  differential)  delays  In  the  received  signals;  the  second 
is  the  development  of  localization  and  tracking  algorithms. 

Previous  work  on  the  multipath  project  included  the  development  of 
differential  delay  estimation  algorithms  as  well  as  theoretical  studies 
resulting  in  bounds  on  the  variance  of  estimating  multipath  delay.  Bounds 
were  also  placed  on  the  accuracy  of  estimating  various  so-called  track 
parameters  —  speed,  depth,  and  range  at  closest  approach  of  the  target  to  the 
sensor  array. 

Under  the  present  contract,  software  was  developed  to  estimate  multipath 
delay  using  real  or  simulated  data  as  Input.  Also  under  the  present  contract, 
tracking  algorithms  were  developed.  More  specifically  the  following  are 
accomplishments  of  this  phase.  ^ 

1.  A  real  data  processing  system  was  developed  for  the  purpose  of 
estimating  mutlpath  delay  as  a  function  of  time.  This  system  was 
composed  of  a  correlogram  generator  and  a  line  tracker.  A 
correlogram  Is  a  picture  whose  rows  consist  of  normalized  correlation 
functions  at  different  time  instances.  The  correlation  functions 
peak  at  the  value  of  the  differential  delay  and  align  to  form  tracks 
In  the  correlogram.  The  line  tracker  pulls  the  differential  delay 
curves  out  of  the  correlogram.  The  correlogram  generator  produces 
SCOT,  ML-,  PHAT-  and  un-normal ized  correlograms.  The  line  tracker 
uses  the  AOEC  line  tracking  algorithm.  This  software  was  developed 
on  the  AP120B  array  processor  on  VAX  11/780  for  high  speed 


computation. 


Two  types  of  track  parameter  estimation  methods  were  developed.  The 
first  makes  use  of  the  measured  value  and  functional  form  of  the 
delay  curve  at  various  points  in  time.  The  other  fits  an  estimated 
delay  curve  to  the  measured  delay  curve,  and  estimates  the  track 
parameters  as  those  giving  the  best  fit  to  the  data.  Target  is 
assumed  to  travel  along  a  constant-depth,  constant-velocity  straight 
line  course,  and  the  parameters  of  the  target's  path  —  depth, 
velocity,  and  closest  point  of  approach  —  are  estimated,  thereby 
localizing  the  target.  The  tracking  algorithms  were  implemented  in 
FORTRAN  77  and  evaluated  using  the  real  data  provided  by  ONR.  The 
results  are  encouraging  in  that  these  methods  allow  tracking  of  real 


underwater  targets. 


TRACK  PARAMETER  EXTRACTION  USING  MULTIPATH  DELAY 
AND  DOPPLER  INFORMATION 

0.  INTRODUCTION 

Acoustic  waves  propagating  In  the  ocean  undergo  various  reflections  and 
refractions  [1].  Because  of  this  multipath  propagation,  acoustic  surveillance 
systems  receive  several  different  delayed  and  doppler  shifted  versions  of  the 
source's  signal.  The  relative  delays  and  doppler  shifts  contain  valuable 
information  which  can  help  localize  and  track  a  target  [2].  Consider,  for 
example,  a  single  multipath  reflection  depicted  in  Figure  0.1.  The  time  delay 
between  the  signals  propagating  in  the  direct  and  reflected  paths  is  a 
function  of  the  target/sensor  range  and  target  and  sensor  depths.  Similarly, 
the  doppler  shifts  of  the  direct  and  multipath  signals  are  functions  of  the 
velocity  and  range  of  the  target.  Therefore,  a  history  of  the  delay  and 
doppler  shifts  can  lead  to  knowledge  of  the  so-called  track  parameters  — 
speed,  range  of  closest  approach,  and  depth. 

In  this  report,  techniques  for  extracting  track  parameters  from  delay  and 
doppler  information  are  presented.  Bounds  are  placed  on  the  accuracy  of  these 
estimates  as  well  as  the  accuracy  of  the  delay  and  doppler  measurements. 

Cases  with  one  target,  one  and  two  sensors  with  and  without  multipath  are 
treated  in  detail. 

The  structure  of  this  report  is  as  follows:  Chapter  one  details  the 
geometry  and  signal  equations  for  a  single  sensor  and  target.  This  chapter 
also  discusses  the  measurement  of  the  delay  and  doppler  shifts  from  the  raw 
data.  Chapter  two  presents  the  delay  and  doppler  formulas  for  a  single  sensor 
and  single  target  in  the  presence  of  a  surface  reflection.  The  accuracy  of 
track  parameter  estimators  using  delay  and  doppler  information  is  also 
discussed  in  Chapter  two.  Chapters  three  and  four  develop  methods  for 
extracting  the  track  parameters  from  delay  and  doppler  information  for  several 
sensor  geometries;  and  chapter  five  contains  some  concluding  remarks. 

There  are  four  appendices:  appendix  one  discusses  the  effect 


Fiqure  0.1  Illustration  of  multipath  due  to  a  surface  bounce  for  the  case 

2  2  1/2 

of  a  single  source  and  sensor.  D  *  (1/C)  [[R  +4zs+4z^z)  -R] 

multipath  delay.  C  is  the  velocity  of  sound  in  water. 


differential  doppler  has  on  the  measurement  of  differential  delay.  Appendix 
two  presents  an  example  of  track  parameter  estimators  applied  to  a  real  data 
set.  Appendix  three  discusses  the  accuracy  of  depth  estimation  using 
multipath  delay.  Finally,  appendix  four  contains  software  listings. 


1.  MULTIPATH  PARAMETERS  FROM  RAH  DATA 

In  this  chapter,  the  extraction  of  delay  and  doppler  information  from  raw 
data  is  discussed.  We  consider  here  the  case  of  a  single  multipath  reflection 
(Fig.  1.1).  The  signal  at  the  sensor  will  be  roughly  given  by 

y(t)  =  x(ad(t).t  -  Dd(t))  +  g(t)x(om(t).tm  -  0  (t))  +  e(t)  (1.1) 

where 

y ( t)  is  the  signal  at  the  sensor 

a  .  .  is  the  doppler  shift  in  the  direct  signal 

dlt) 

a  (t)  is  the  doppler  shift  in  the  multipath  signal 
m 

D  .(t)  is  the  delay  in  the  direct  signal 
d 

D  (t)  is  the  delay  in  the  multipath  signal 
m 

g(t)  is  the  gain  of  the  reflected  signal 
e(t)  is  the  sensor  noise 

x(t)  is  the  target  signal 

It  is  assumed  that  e  is  uncorrelated  with  x;  the  doppler  shifts  are  small: 

“d*  “m  ”  1  »  and  the  delays,  doppler  shifts  and  gains  are  slowly  changing 

with  time:  ||!jr| ,  j||| ,  |||j  «  1;  and  0  «  |g|  <  1. 

1.1  DELAY 

Assuming  a.  =  a  =  1  ,  the  Cramer-Rao  lower  bound  on  the  variance  of  the 
d  m 

minimum  variance  estimator  of  the  differential  delay,  D  =  ID  .  -  D  I  ,  is 

d  m 

reasonably  small  for  signals  with  a  large  bandwidth-delay  product  [3].  When  x 
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Figure  1.1  The  received  sensor  signal  in  the  presence  of  a  multipath 
refl ection. 


and  e  are  uncorrelated  lowpass  white  noise  processes  with  bandwidth  w,  the 
variance  in  estimating  D  from  the  signal  y(t)  =  x(t)  +  g.x(t-D)  +  e(t)  is 
bounded  by. 


Var(D)  >  JZ-rL— 

,  SNR  «  1,  WD  »  1 

(1.2a) 

TwV  SNR^ 

VAR(D)  >-j^l-g2) 

,  SNR  »  1,  WD  »  1 

(1.2b) 

w  g 

where  the  signal-to-noise  ratio  (SNR)  is  the  ratio  of  the  power  contained  in 
the  signal  and  noise.  The  maximum  likelihood  (ML)  estimator  achieves  this 
bound,  but  involves  a  computationally  intensive  minimization  of  a  complicated 
cost  function  [2]. 

Maximizing  a  simplified  version  of  the  ML  cost  function  results  in  the 
autocorrelation  method  of  multipath  delay  estimation:  the  delay  is  chosen  as 
the  value  of  i  which  maximizes  | Ryy ( t ) I  ,  the  absolute  value  of  the 
autocorrelation  function  [2].  The  autocorrelation  function  is  relatively  easy 
to  compute,  and  Is  commonly  used  for  delay  estimation  [4], 

The  variance  and  bias  in  estimating  delay  using  the  autocorrelation 
method  based  on  data  over  a  time  period,  T  are  given  in  [2]  by 


Var,D*,  =  -  ^sM)2 


SNR  «  1,  HD  »  1 


(1.3a) 


BIAS  =  < 


2 

w  u 


gs. 


-r 


(1.3b) 


Note  that  for  low  SNR,  the  ML  estimator  and  the  autocorrelation  estimator  have 
essentially  the  same  variance.  Therefore,  for  sufficiently  large  bandwidth, 
low  SNR  signals,  as  is  usually  the  case  in  underwater  acoustic  surveillance, 
the  correlation  estimator  will  perform  almost  optimally. 


Tracking  Differential  Delay  D(t) 

For  Oft)  «  1.  and  a  =  3  1.  the  autocorrelation  function  and  power 

»  m  a 

soectrum  of  y(t)  are  given  by 


Rvv(x)  =  (1+g  )RYY(x)  +  gRYY(x-D)  +  gRYY(T+0)  +  Rpp(x) 


S  yy  ( cjj )  =  (1  +  g  +  2gcoscaD)SYY  (  oj  )  +  Sp-(uj) 


(1.4a) 


(1.4b) 


From  (1.4),  it  is  easily  seen  that  when  the  width  of  Rxx(x)  <<  1  ,  peaks  of 

Ryy{t)  will  occur  at  the  value  of  the  multipath  delay  as  well  as  at  zero 

del  ay . 

When  the  SNR  is  low,  and  D(t)  is  slowly  changing  with  time,  D(t)  is  most 

often  found  using  a  correlogram  —  a  2-D  gray  level  image  whose  i-th  row  is 

the  correlation  function  taken  at  t. ,  E(Y(t.)Y(t.+t  ))  See  Figure  A3.1  for 
example.  Peaks  in  the  correlation  functions  comprising  a  correlogram  will 

merge  to  form  lines  tracing  the  history  of  D(t). 

Doppler  Effects 

Since  the  direct  and  multipath  signals  are  arriving  from  different 
directions,  a  moving  target  will  cause  a  relative  doppler  shift  between  the 
two  signals.  The  autocorrelation  of  the  sensor  signal  will  become 

Ryy  ( t ,  x )  =  R^Ht-xla^]  +  g^R^^[(t-T)aml 

+  9RxxCod t’  V  "  01  +  ^XX^d*’  V  +  D]  +  Ree(t'T)]  (1,5) 


where 


Rxx[a,b]  =  E[X(a)X(b)]  ,  and  r^t]  =  E[X(t)X(t+T)] 

As  shown  in  Appendix  1,  the  secondary  peaks  in  Ryy(x)  which  are  used  to 

estimate  the  delay  are  reduced  as  the  differential  doppler,  a .  -  a  , 

dm 

increases.  That  is,  the  direct  and  multipath  signals  become  decorrelated  as  a 


sww 


»  "  >  “  v  -  W  * 


result  of  target  motion. 


The  differential  doppler  is  greatest  when  the  difference  in  the 
directions  of  arrival  of  the  direct  and  multipath  signals  is  greatest.  This 
occurs  at  CPA,  the  point  of  closest  approach  of  the  target  to  the  sensor 
array.  The  fading  of  the  correlogram  lines  near  CPA  in  Figure  A2.1  is  in  part 
due  to  this  differential  doppler  ef'ect. 

1.2  DOPPLER 

Assuming  D_,  =  D  =  0  ,  substituting  t'  =  log  t  into  (1.1)  gives 
d  m 

Y(t')  =  X(t'+loga.)  +  X ( t '  +1  ogct  )  +  e(t')  (1.6) 

a  m 

Note  that  since  e(t)  is  white  gaussian  noise,  so  is  e(t').  Therefore, 
estimating  differential  doppler  in  the  absence  of  delay  is  equivalent  to 
estimating  differential  delay  in  the  absence  of  doppler. 

The  ML  estimator  of  doppler  shift  achieves  the  C-R  bound  of 

Var(log^-)  =  — _L_  ,  snr  «  1  (1.7a) 

“m  W  V  SNR4 

Var(log^)  =  -^-(l-g2) ,  SNR  »  1  (1.7b) 

“m  W,3g2 

where  W  is  the  bandwidth  of  the  log  sampled  signal,  X ( t ’ ) .  The 
autocorrelation  method  of  determining  differential  delay  estimation,  looking 
for  the  peak  of  E ( Y ( t ' )Y(t'+loga))  ,  i.e.,  the  peak  of  E(Y\t)Y(at) )  ,  has 
variance 

Var(loga)  =  -34-y  -U-  ,  SNR  «  1  (1.8) 

W*  Jgg  SNR4 

where  it  is  assumed  that  loga  •  W  <<  1  . 

Tracking  g(t) 


The  doppler  ratio,  a  =  — —  is  often  seen  using  dopplergram.  Analogous 

ad 


9 


to  a  correlogram,  a  dopplergram  is  a  2-D  image  whose  i-th  row  is 
E (Y ( )Y (ctt^ ) )  .  Lines  in  the  dopplergram  show  the  time  history  of  a(t)  . 

Delay  Effects 

As  differential  doppler  decorrelated  the  multipath  and  direct  linearly 
sampled  signals,  differential  delay  will  decorrelate  the  multipath  and  direct 
log  sampled  signals: 

E[Y(t)Y(a*t-D)]  “  gRxx(0)  <  gRxx(0)  (1.9) 

* 

where  a  is  the  correct  value  of  the  doppler  ratio.  Differential  delay  is 
greatest  at  CPA.  Consequently,  one  would  expect  fading  of  differential 
doppler  lines  around  CPA. 

Absolute  Doppler 

As  a  final  note,  absolute  doppler  information  can  be  obtained  by 
following  the  motion  of  narrow  band  components  in  the  spectrum  of  the  received 
signal.  A  2-D  image  made  of  spectra  taken  at  different  times,  a  spectrogram, 
is  often  used  for  this  purpose. 


2.  DELAY  AND  DOPPLER  EQUATIONS,  AND  ACCURACY  ANALYSIS 

This  chapter  develops  basic  delay  and  doppler  formulas  used  in  finding 
track  parameters  for  measurements  of  differential  delay  and  doppler;  bounds 
are  placed  on  the  variance  of  track  parameter  estimators  using  these 
measurements.  It  is  assumed  here  that  the  ocean  is  a  homogeneous  medium  with 
a  reflecting  boundary  at  the  surface  (z  *  0). 

2.1  DELAY  AND  DOPPLER  EQUATIONS 


With  a  sensor  at  (Xsl  Ys>  Z$)  and  a  target  at  (vt,  Y^,  Lj.)  ,  the 
direct,  multipath,  and  differential  delay  between  direct  and  multipath  signal 
are  given  by  (see  Figs.  1.1,  and  2.1): 


Dd(t)  =  [(Xs-Vt)2  +  (Ys-Yt)2  +  (Z$-Zt)2]1/2/C 

(2.1a) 

Dm(t)  =  C(xs-Vt)2  +  (YS-YT)2  +  (Zs+Zt)2]1/2/C 

(2.1b) 

0(t)  £  Dm(t)  -  Dd(t) 

(2.1c) 

where  C  is  the  speed  of  sound  in  water.  The  differential 
differential  delay  curvature  are: 

delay  rate  and 

iH  =  (y )  r! _ !_i 

3t  r2  lVX  Xs  ‘■D_  D.J 

L.  m  d 

(2.2a) 

32D  V2  1  1  (Vt-X  )  .  . 

-r  *  rArr  -  rrl  *  — 3 

(2.2b) 

m 


The  doppler  shift  of  a  source  frequency,  f$  is  given  by 

f  =  V1  +  h  --  V6 


(2.3) 


where  f  is  the  received  frequency,  and  Vr  is  the  velocity  of  the  target  along 
a  line  connecting  the  target  and  sensor,  and  a  is  the  doppler  shift.  The 
direct,  multipath,  and  differential  doppler  shift  between  the  direct  and 
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TOP  VIEW 


Figure  2.1 


View  from  abovo  the  ocean  of  the  target  and  sensor.  The  target 
is  located  at  (X^-Vt,  Y^. ,  Z^)  and  the  sensor  is  located  at. 


12 


multipath  signals  are  given  by 


V{V t-X_)  3D j( t) 


V(Vt-Xj  3D  ( t) 

6m(t)  s  1  +  'TZT1'  =  1  +  ~fr- 


C  0 


m 


*<«  *  <■<*>  -  «dll)  ■  ^r1 


(2.4a) 


(2.4b) 


(2.4c) 


The  differential  delay  between  direct  paths  from  the  target  to  two 
sensors  is  given  by 


D 


12 


(2.5a) 


where  Drfl  and  Dd2  are  the  delays  between  the  target  and  the  individual 
sensors.  The  differential  doppler  between  the  direct  paths  of  the  two  sensors 
is  again  given  by  the  differential  delay  rate: 


aP12 

*12  "  <Sdl“'Sd2  =  IT 


(2.5a) 


(2.5b) 


2.2  ACCURACY  ANALYSIS 

It  is  of  interest  to  know  how  well  the  track  parameters  can  be 
estimated  from  multipath  delay  Information.  The  C-R  bound  can  be  used  in 
conjunction  with  the  above  formulas  to  put  a  lower  bound  on  the  variance  of 
any  unbiased  estimators  of  the  track  parameters. 

C-R  BOUND  CALCULATION 


It  Is  assumed  that  data  is  available  during  the  time  period  (-T,T).  Over 
this  time  period,  the  data  are  nonstationary.  Since  the  C-R  bound  can  be  used 


only  with  stationay  data,  for  purposes  of  determining  a  bound  on  the  variance 
of  the  track  parameter  estimates,  it  is  assumed  that  estimates  of  track 
parameters  are  made  from  sub- interval s  of  time  of  length  At  (over  which  the 
data  are  assumed  stationary)  and  linearly  combined  to  form  a  single  estimate 
of  the  track  parameters.  In  this  case  the  variance  in  estimating  a  track 
parameter,  e  is  given  by. 


A 

and  0  is  given  by, 

A 
8 

i=-N  <jg(i)'  i*-N 


2  =  r  t  i-  ■ 

*9  1=-N  <j2(  i ) 


-1 


(2.6a) 


N 

(  t 


a  ( i )  ^ 

4lL)/(  y 


«?nn 


(2.6b) 


where  the  variances,  <j2(i)  are  bounded  below  by  the  C-R  bound, 

0 

°8{1)  >  TTTT  (2*7) 

9 

where,  1(1)  is  the  Fisher  information  matrix  element  corresponding  to  the 
9 

parameter  a  at  time  t^  . 

The  variance  in  estimating  e  from  time  delay  measurements  over  the 
interval  (-T,T)  is  therefore  bounded  by, 

aQ  >  -  (2.8) 

l  I  (i) 

1--M  9 

which  can  be  approximated  by, 

a2  >  [  JT  I9(t)dt]_1  (2.9) 

-T 

for  large  N.  The  Fisher  information  matrix  element,  I  (i)  ,  can  be  evaluated 
,  8 
exactly  using  methods  described  In  Appendix  3.  These  expressions  are 

algebraically  complicated,  and  difficult  to  manipulate.  If  it  is  assumed 

that  all  track  parameters  but  8  are  known,  then  I  (t)  is  given  simply  by, 

9 


where  I Q ( t; )  Is  the  Fisher  Information  matrix  element  for  the  delay,  evaluated 
in  section  1.  Note  that  the  Fisher  information  above  is  larger  than  that 
calculated  assuming  the  other  track  parameters  are  unknown.  Therefore,  the 
resulting  bounds  on  the  variance  of  the  track  parameter  estimates  will  not  be 
tight. 

Using  (2.1),  (2.10),  and  (2.9),  a  lower  bound  can  be  placed  on  the 
variance  of  track  parameter  estimates. 


’» >  l  Hit1)2 


a  *■  J  ^  39  >  D 

-T 


(2.21) 


where,  from  (1.2a), 


ID(t)  -  w3g2SNR2/3 


for  low  SNR  broadband  signals. 


DEPTH  ESTIMATE  VARIANCE 


Assuming  other  track  parameters  are  known,  the  Fisher  information  of 
the  depth  estimate  based  on  delay  measurements  is  given  by  (2.10)  as. 


(Z,-Zc)2  (Z,+Zc)2 
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Using  (2.21)  and  integrating, 

1  w3 g2 SNR2  (ZT-Z<.)2  ,  (ZT+Z«.)2  i 

T“  <  (-T— )  *  ^n"1  tan"1  ~ 


(2.13) 


,2  _2  t=X_+VT 

+  2  lllll  tan-l  I_} 

°d  ad  t=Xt-VT 


where 


(Y2  +  (2t-Zs)2)1/2 
(Y2  +  (ZT+ZS)2)172 
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In  the  case  of. 


X-  t  VT  <<  a.  .,  a  ,  a 
r  a  m  + 


2  2  2 
YT  »  ZS*  A 


(2.14) 


(2.15) 


reduces  to. 


1  w3q2SNR2  4Z2 

~r~  *  •  fcr)  •  2T 

a7  oC  Y  T 

LJ  1 


(2.16) 


Notice  that  this  expression  is  proportional  to  T,  i.e.  the  more  data 
available,  the  lower  the  variance.  Also  notice,  as  the  depth  of  the  sensor 
increases  relative  to  the  range  at  CPA,  the  minimum  variance  is  decreased. 


When  (2.15)  and. 


1 X±VT |  >>  ad,  o^,  a+ 


(2.17) 


hold,  i.e.  most  all  of  the  delay  curve  is  available, 

i  wVcup2  *Z2 

4-  *  (  --H— HvHr1)  (2-18) 

a7  3C  T 

i 

In  this  case,  the  variance  in  estimating  the  target  depth  is  dependent  on  the 
sensor  depth,  the  deeper  the  sensor,  the  better  the  depth  estimate;  and  the 
target  velocity  and  CPA  range,  closer  slower  moving  targets  give  better  depth 
estimates. 
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Y  OFFSET 


With  (2.11)  and  (2.1), 
known,  given  by. 


in  the  case  of  other  track  parameters 


(2.19) 
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a+  a+  't=Xt-VT 


If  only  data  near  CPA  Is  available,  i.e.  (2.14)  applies,  then  under  (2.15), 
—  reduces  to. 
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(2.20) 
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In  other  words,  the  ay  Is  reduced  for  increasing  Z$,  Zy,  and  T  and 
decreasing  target  CPA  range. 

When  a  lot  of  data  are  available,  i.e.  when  (2.17)  is  applicable, 
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av  becomes , 
tT 
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(2.21) 


In  this  case,  more  data  no  longer  significantly  improves  the  estimate.  Note 
here  that  the  slower  the  target  moves,  the  broader  the  delay  curve,  and  the 
better  the  estimate. 

X  OFFSET 


Using  (2.10)  and  (2.11)  the  variance  of  the  X  offset  can  be  bounded  by. 


w3q2SNR2, 


Vvt 


-°dtan'1  77  _  amtan  1  ^~+  2a+tan  1  J"l|T-xL 
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(2.22) 


Under  the  conditions  (2.14)  and  (2.15), 

2  2 

1  .  wV$NR2  4  „2,3,ZsZT,  , 
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Note  here  the  heavy  dependence  of  on  V,T,  and  YT  :  the  more  data 

XT  1 


(2.23) 


available  near  CPA,  and  the  more  peakld  the  delay  curve  is,  the  easier  it  is 
2 

to  estimate  <,*  • 

*T 

When  (2.27)  applies. 
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(2.24) 


In  this  case,  the  estimate  variance  is  no  longer  as  dependent  on  the  shape  of 
the  delay  curve  or  the  amount  of  data  available. 


VELOCITY 


Again,  using  (2.10)  and  (2.11),  a  bound  on  the  variance  of  the  velocity 
estimate  can  be  found  as. 


+  2XT(r  log(°d  *  T  )  +  T  1og(am  +  c  )  "  °+  1og(a+  +  T  ' 

p  i  i  i  t=Xt+VT 

+  X^-a^tan'1  - —  a_  tan'1  —  +  2  tan"  — )  1  ' 

d  ad  m  am  a+  °+  t=Xt-VT 

When  (2.14)  and  (2.15)  are  satisfied. 


1  .  w3g2_SNR^  .(-2V5T5+XTV4T4+  £  XyV3T3)  2(^^) 


LEYVA'S 


Again,  the  estimate  variance  is  sensitive  to  V,T,Y^.  ;  lower  variance 
estimates  are  obtained  from  data  near  CPA  where  the  delay  curve  is  more 
peaked. 

If  instead  of  (2.14),  (2.15)  is  applicable,  then 


2*w3g2SNR2 


(Zs  Zl) 


(2.27) 


l 


Notice,  in  this  case,  the  variance  of  the  velocity  estimate  increases  with  the 
"sharpness"  of  the  delay,  and  not  the  "flatness"  as  was  the  case  in  (2.26). 

SUMMARY 

In  summary,  provided  the  delay  estimates  are  good,  i.e.  the  banwidth  of 
the  signal  is  large  enough,  and  provided  the  aperture  size  of  the  array  is 
large,  i.e.  Z$  is  large  enough  compared  to  the  CPA  range  of  the  target, 
reasonably  small  lower  bounds  on  the  variance  of  the  track  parameter 
estimates  are  obtained. 

Tighter  bounds  on  the  track  parameter  estimates  as  well  as  bounds  for 
estimates  based  on  inter-sensor  and  multipath  delay  information  from  two 
sensor  arrays  are  presented  in  appendix  3. 
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3.  POINT  METHOD  OF  TRACK  PARAMETER  EXTRACTION 


Techniques  for  extracting  track  parameters  from  correlograms, 
dopplergrams,  and  spectrograms  are  presented  in  this  chapter.  The  parameters 
of  interest  are  V,  the  velocity  of  the  target;  zj,  the  depth  of  the  target; 
Rcpa»  the  radius  of  closest  approach  to  the  sensor  array;  and  9  ,  the  bearing 
angle  relative  to  the  array. 


3.1  SINGLE  SENSOR 


The  case  of  a  single  sensor  at  (0,0,  Z$)  and  a  single  target  at 
(Vt.Yy,  Z.p)  with  one  surface  reflection  is  considered  first.  Figure  3.1 
shows  the  correlogram  for  this  case.  Differential  delay  rate,  estimated  from 
the  dopplergram  is  plotted  in  Figure  3.2.  These  measurements  lead  to 
estimates  relating  R^,  V,  and  Zj  .  Measurements  of  the  doppler  shift  of 
spectral  lines  give  velocity  as  well  as  estimates. 


Delay  Measurements 


When  the  target  is  far  away  from  the  sensor,  the  differential  delay 
becomes  less  dependent  on  the  relative  y-axis  positions  of  the  sensor  and 
target.  As  |t|  +  •  ,  a  relationship  between  Zj  and  V  can  be  developed. 


From  (2.1)  for  |t|  »  1, 
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Since  for  |t|  »  1,  (Vt)2  »  'll  . 


ffi 


Performing  a  least  squares  fit  of  p/t  to  the  tails  of  the  measured  D(t) 
yields  the  rel ationship: 


-  C 

T~  =  p  7T~ 
s 


(3.2) 


where  p  =  min[0(t)  -  -^]  for  t  >>  1 


A  relationship  between  Rro.  and  ZT  can  be  found  by  examining  0(0),  a  readily 

2  2  1/2 

available  quantity.  From  (2.1)  and,  RCpA  =  Dd(0)  =  [YT+(ZS  -  Zy)  ] 


COO)  -  [>& A  *  «sZt]1/2  -  “CPA 


(3.3) 


Solving  (3.3)  for  R 


CPA  ’ 

n  _  2ZSZT  CD(0) 
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(3.4) 


Using  a  velocity  estimated  by  other  means,  Rcpa  and  Zj  can  be  obtained 
using  (3.2)  and  (3.4)  with  a  few  measurements  of  D(t). 

Doppler  Measurements 


By  following  doppler  shifts  of  lines  in  the  spectrogram,  V  and  R^pA  can 
both  be  estimated.  Figure  3.3  shows  the  spectrogram  of  a  constant  velocity 
target  emitting  a  single  strong  spectral  line.  Using  (2.3)  and  (2.4), 


f  •  V1+  cr)  =  y 1+ 


V(Vt-Xs) 


(3.5) 


as  t  +  D  ,  D 


d  m  C 


and 


f  >  f_ 


t1  *  bf, 


(3.6) 


Similarly,  as  t  +  +»  , 
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Figure  3.3  Doppler  shift  vs.  time  for  a  constant  velocity  source 
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from  (3.6)  and  (3.7)  an  expression  for  the  velocity  is  obtained, 


fs  =  2  (fmax 


^min) 


V  =  C 


f  -  f  . 
max  min 


max  min 


af 

Range  information  can  be  found  by  examining  at  CPA: 

3t 


fs  3Vr| 


CPA  C  3t  CPA 


(3.10 


from  (3.7) 


Hi  s  f  I  =  li  vi_ 

atjCPA  s  3t2  lCpA  TT  rtCf>A 
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The  time  of  maximum  differential  doppler  is  easily  obtained  from  the 
dopplergram,  and  relates  the  direct  and  multipath  CPA  ranges.  The  maximum 
differential  doppler  time,  t*.  is  given  by  (see  (2.4)), 


D(t)  =  Dm(t*)  -  Dd(t*)  =  0 


(3.12 


Using  (3.2) 
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Therefore, 


**  =  ^IDm(t*)Dd(t*) 


(3.13) 


5(t*)  =  D(t*) 
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(3.14) 


Using  (2.1),  (2.2)  and 


«d  ■  V  ♦  <W 


",  -  A  *  <  w2  • 


the  time  of  maximum  differential  delay  is  given  by. 


"  W  [Rd  +  4  +  <*d  +  +  32R2dR2J1/2]1/2 


(3.15) 


Summary 


To  summarize,  from  the  correlogram  and  spectrogram,  all  track  parameters 


can  be  estimated. 
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p  =  Mi n[D(t)  -  p/t]' 


(3.17a) 


ZT  =  p  v  *  7T 
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(3.17b) 


(3.17c) 


£  i  CD(O) 

KCPA  ‘  COWT  '  1 

where  f  .  f  .  .  »  are  measured  from  the  spectrogram;  l  and  D(0) 

max’  min’  at|t=RCPA  p 

are  measured  from  the  correlogram. 


3.2  TWO  SENSOR  TRACK  PARAMETER  ESTIMATION 


With  the  addition  of  a  second  sensor,  the  track  parameter  estimation 
capability  of  an  array  is  greatly  improved.  First,  V,ZT  and  Rcpft  estimates 
from  each  sensor  can  be  combined  to  form  estimates  of  lower  variance.  Second, 


measurements  of  intersensor  delay  and  doppler  lead  to  further  estimates  of 
V,  Zy,  RCpA  ,  and  possibly  0  -  the  bearing  angle  relative  to  the  array. 


Two  commonly  used  two-sensor  arrays  are  considered  here,  the  vertical 
array  and  the  horizontal  array.  The  arrays  are  shown  in  Figure  3.4; 
theoretical  correlograms  in  the  presence  of  multipath  are  shown  in  figures  3.5 
and  3.7. 


The  Vertical  Array 


The  vertical  array  depicted  in  Figure  3.4  has  two  sensors  placed  one 
above  the  other.  This  array  is  radially  symmetric,  and  therefore  can  only 
estimate  Zy,  V,  and  RcpA  .  Separate  estimates  of  these  parameters  can  be 
made  from  multipath  information  at  each  sensor  as  well  as  inter-sensor  delay 
and  doppler  shift  data.  Cross  correlation  and  cross  differential  doppler  are 
shown  in  Figure  3.5. 


As  can  be  seen  in  Figure  3.6,  the  delay  between  the  direct  and  multipath 
propagation  paths  in  the  single  sensor  case  is  equivalent  to  the  Intersensor 
direct  path  delay  in  the  2-sensor  vertical  array.  If  the  sensors  of  the 
vertical  array  are  placed  at 


SI:  (0,0, Zsl) 
S2:  (0,0, Zs2) 
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Figure  3.4  Vertical  and  horizontal  2-sensor  arrays 
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3.6  The  equivalence  of  a  two  sensor  vertical  array  in  the  absence  of 
multipath  and  a  single  sensor  array  in  the  presence  of  multipath 
can  be  seen  in  the  above  diagram. 


where 


's2 


>  Z 


si 


then,  making  the  substitution 


Z 


S 


Zs2-Zsl 


(3.18) 


in  (3.17)  will  result  in  track  parameter  estimates  from  inter-sensor  delay  and 
doppler  estimates.  Note  that  Rqpa  in  this  case  is  measured  as  the  distance  of 
closest  approach  to  the  deepest  sensor.  Also  note  that  Zy  is  the  depth 
relative  to  the  midpoint  between  the  sensors. 


The  various  estimates  of  V,  ZT  and  R„.  (sensor  1,  sensor  2,  and  inter- 

T  CPA 

sensor  measurements)  should  be  combined  so  as  to  minimize  the  variance  of  the 
final  estimates.  The  method  described  in  (2.6)  is  optimal  over  all  linear 
combinations  of  the  different  estimates.  In  this  case,  more  weight  is  given 
to  the  smaller  variance  estimates.  If  the  multipath  reflection  is  weak,  the 
inter-sensor  measurements  are  more  important.  Otherwise,  the  estimate  made 
from  the  array  with  the  larger  aperture  size  will  have  the  least  variance  and 
therefore  more  weight.  (The  aperture  size  is  2Z.  for  the  ith  sensor,  and 
Zc  -Zc  for  inter-sensor  measurements). 

i 

The  Horizontal  Array 

The  horizontal  array,  depicted  In  Figure  3.4,  consists  of  two  sensors  at 

the  same  depth.  As  this  array  is  not  radially  symmetric.  Intersensor 

multipath  information  leads  to  estimates  of  bearing  angle  as  well  as  and 

V.  These  estimates  can  be  combined  with  estimates  from  data  at  Individual 

sensors  to  form  estimates  of  Zj  ,  the  depth  of  the  target;  Rcpft  ,  the  radius 

of  closest  approach;  V,  the  target  velocity;  X  and  Y  ,  the  X  and  Y  axis 

c  c 

crossings  of  the  target;  and  9  ,  the  bearing  angle.  Due  to  symmetries  in  the 
array,  however,  9  and  Yj.  can  be  only  determined  to  within  a  sign. 


The  horizontal  array  has  sensors  located  at 
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The  target  is  moving  by  the  sensors  at  angle  0  relative  to  the  X  axis,  at  a 
velocity  V,  and  a  depth  of  . 

The  position  of  the  target  is  given  by. 

Target:  (Vtcose  +  X^.,  Vtsine  +  Y^,  Z^). 

So  the  target  is  at  CPA  when  t=0,  define,  ctne 

Yt  =  X  etna  . 

Note  that  a  target  travelling  along  this  trajectory  will  have  axis  crossings 
of 

Xc  =  XT(1  +  ctn2e) 

Yc  =  YT(1  +  Tan^a) 

From  (2.5)  the  inter- sensor  delay  is  given  by 

2  2  2  1/2 
D12( t)  =  [(Vtcose  +  XT  +  x$)  +  (Vtsine  +  YT-YS)  -  (Zs-ZT)£]A,e 

-[(Vtcose  +  XT  -  x$)2  +  (vtsine  +  YT-Y$)2  +  (Z$-ZT)2]1/2  (3.19) 

Figure  3.7  shows  the  cross- correlogram  and  differential  doppler  for  this  case. 

By  examining  the  far  range  behavior  of  (3.19),  an  estimate  of  9  can  be 
obtained.  When  the  target  is  far  from  the  array,  signals  emitted  from  the 
target  arrive  at  the  array  from  essentially  the  same  direction,  0  —  the 
bearing  angle.  The  intersensor  time  delay  will  then  be 

2X 

0  ( 1 1 1  •)  -±  cose  , 

The  bearing  angle,  a  can  therefore  be  found  as: 
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Examining  the  intersensor  delay  at  CPA,  (at  t=0),  yields  an  expression 
relating  and  R^pA  .  Recall, 


'CPA 


(Xj  +  +  (ZT  -  Zsn 


2,1/2 


(3.21) 


Using  (3.17) 


,1/2 


0012(0)  =  CRCPA  +  ^RCPA  “  2XTXS] 
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(3.22) 


Sovling  for  R 


CPA  * 


rCpa  *  C(C2D22(0)  -  2XtXs)(C2D22(0)  -  6XtXs)]1/4 


(3.23) 


When  the  relative  sensor  delay  Is  zero,  the  target  is  between  the 
sensors,  at  (0,  Y  ,  Z^.)  .  Examining  the  relative  delay  rate  at  this  point 


yields  an  expression  for  Y  .  Define 

V* 


ty  &  t:012(t)  =  0 
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(3.24) 


Recall , 
3D 
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C  0,  CD, 


+  ^fj^Vtsine  +  VT-XS)  -  —^(Vtcose  +  XT  -  X$) 


(3.25) 


where  D1  and  02  are  the  delays  from  the  target  to  sensor  1  and  sensor  2. 
When  D12  *  0,  0^  =  D2  ,  and  the  target  Is  at  (0,Yc>  Lj.)  .  Therefore, 
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At  t. 


c2°i  *  (ZT  -  Zs)2  +  Y 


and,  solving  for  Y 


Yc  *  [X|[(2-~C~~— )2-  1]  -  (Zt-Zs>2]1/2 

C0(t  ) 

Yc 


(3.27) 


When  the  differential  delay  between  sensors  is  at  a  maximum,  the  target 
is  at  (X^ ,  0,  Lj.)  .  Looking  at  D^(t^  )  .  where 


^  Max  1 012 ( t)  I  , 

yields  an  expression  for  X  .  At  t„  ,  0,_(t)  is  given  by 

c  Xc  12 

CD12(tX2)  =  [{Xs+Xc)2  +  (ZS-ZT^1/2  -  ^xs“xc)2  +  {ZS_ZT)]1/2 


(3.28) 


(3.29) 


Solving  for  X  , 
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cV-U  ] 

- 5 - —  +  (Xs+*z  s^V^^X  >) 


4Xs  -  C 

The  sign  of  X  can  be  determined  by  looking  at  the  times  of  CPA  of  the 
c 

individual  sensors;  X  will  have  the  sign  of  the  sensor  location  that  had  the 

c 

first  CPA.  Note,  however,  if  Z  "  ZT  or  the  bearing  angle  is  close  to  0, 

2X  s 

then  and  xc  cannot  be  determined  accurately  in  this  manner. 


Velocity  can  be  determined  by  measuring  various  times  of  CPA  and  axis 
crossings.  Four  position- time  pairs  can  be  fairly  easily  measured. 


(X  ,0,Z^  at  t  s  max  D  (t) 
c  T  xc  t  12 

(0,  Y  ,  Z_)  at  t  =  t:D  (t)  =  0 
c  T  yc  12 


MW 
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(a+ZXesinecos9,0  +  Z¥  cos  9,ZT)  at  t,  =  max  D  it) 
s  *s  1  kCPA  2  t  1 

where  D^t)  and  02(t)  are  the  differential  delays  between  the  direct  and 
multipath  signals  at  sensor  1  and  sensor  2;  and  a,  s  are  constants.  Two  of 
the  many  estimates  of  velocity  are. 


v  -  cx|  -  ^)1/2 


(3.31a) 


(3.31b) 


To  summarize,  estimates  of  9,  Y  ,  Xc,  R^pA  ,  and  V  can  be  obtained  as 
functions  of  the  inter-sensor  delay  and  doppler  (delay  rate),  Z and  sensor 
coordinates.  These  esitmates,  in  combination  with  single  sensor  multipath 
data  completely  describe  the  track  parameters  of  a  target  moving  along  a 
straight  line  past  the  array. 


t  V 


4.  PARAMETRIC  FIT  fCTHOOS  OF  TRACK  PARAMETER  ESTIMATION 


The  parameter  extraction  techniques  presented  in  the  previous  chapter 
relied  on  knowledge  of  D( t)  and  derivatives  at  specific  points  in  time, 
namely  | t | <<1  (near  CPA)  and  |t|>>!  (far  from  CPA).  Often  times,  D ( t )  is  not 
known  accurately  in  these  regions  of  time  (see  appendix  1  or  section  1.1), 
and  the  parameter  extraction  techniques  described  in  the  previous  chapter  can 
not  be  used.  This  chapter  develops  parameter  extraction  techniques  based  on 
all  available  delay  information,  not  just  that  near  and  far  from  CPA.  Two 
different  techniques  are  developed  and  applied  to  single  sensor,  two  sensor 
vertical,  and  two  sensor  horizontal  arrays. 

As  the  functional  form  of  the  delay  is  known  (see  section  2),  a 
parametric  fit  of  a  model  D  can  be  made  to  the  measured  D.  The  parameter 
estimates  would  then  be  chosen  as, 


(4.1) 


(Xt,Yt,  Zt,  V)  =  min  d(D,D) 

(xt,yt,zt,v) 

where,  CD  is  the  measured  delay,  and  recall 

CD  =  [ (XT-Vt)2  +  Yj  +  (Zt-Zs)2]1/2-  [ (XT-Vt) 2  +  Y2  +  (ZT+Z$)2] 1/2 


for  some  distance  measure  d. 

The  minimization  (4.1)  is  over  a  cost  function  which  is  non-convex  in 
the  track  parameters  for  common  distance  measures.  Therefore,  exhaustive 
search  methods  must  be  used  in  finding  a  solution.  If  each  of  the  four 
parameters  were  quantized  to  be  one  of  n  values,  and  there  were  m  data  points 
available,  the  minimization  (4.1)  would  require  0(mn**4)  computations;  about 
25  hours  compute  time  for  m=500,  n=100,  on  a  1  mflop  computer.  To  overcome 
this  computational  burden,  (4.1)  can  be  reformulated  so  that  functions  of 
(Xj,Y^.,Zy,V)  appear  as  linear  coefficeints  in  a  least  squares  minimization, 
requiring  only  0(4m)  computations,  and  less  than  one  second  compute  time  on 
the  same  computer.  Two  linear  least  squares  formulations  of  (4.1)  are 
developed  below.  The  case  of  a  single  sensor  listening  to  a  target  in  the 
presence  of  a  multipath  reflection  is  considered  first,  and  later  generalized 
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to  the  case  of  two  sensors 
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4.1  EQUATION  ERROR  APPROACH 


This  method  is  an  equation  error-like  approach.  Recall, 


o  »  P-Q, 


(4.2) 


where. 


P  =  [(Xrvt).2  +  Yt2  +  (Zt+Zs)2]1/2/C 
Q  =  [(XT-Vt)2  +  Yt2  +  (Zt-Zs)2]1/2/C 

What  is  desired  is  a  relation  in  terms  of  0,  p  ,  and  Q  so  that  functions  of 
track  parameters  appear  as  linear  coefficients  which  can  be  fit  using  least 
squares  techniques.  Manipulating  (4.2)  yeilds  the  following  relationships. 


D2  3  P2  +  Q2  -  2PQ 


(D2  -  P2  -  Q2)2  =  4P2Q2 


(CD)4  -  4V2(cDt)2  +  8VXT(C2D2t)  -  4(YT2+ZT2+Z2+XT2) (CD)2 


(4.3) 


(4.4) 


(4.5) 


2  2 

+  16Zj  Zj  =  0 


With  values  of  D(t)  specified  at  m  time  instances,  solving  for  the 
coefficeints  of  (4.5)  which  minimize. 


P  =  min  i f-fi 2  ,  f  =  $A 
A 


(CD(t1)t1)2  C2D(t1)2t1  C20(t1)2  1 


(4.6) 


(CD(tn,)tJ2  C20(tm)2tm  C2D(tm)2  1 
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is  a  standard  linear  least  squares  problem. 


A  = 

al1 

a2 

f  = 

(CD(t1) )4 

• 

• 

• 

a3 

a4 

• 

(CD(tJ)4 

m 

The  optimal  set  of  coefficeints,  A*  is  given  by, 

A*  =  UT*fl  *Tf  (4.7) 

and  the  resulting  track  parameter  estimates  are. 


r  \rh 
-  7  *2//3l 


1  , 

T  (a3 


VZ. 


A" 

4T“ 

s 


i_  a/2 
16  J 


SENSITIVITY  ISSUES 


(4.8) 


Upon  closer  examination  of  (4.1)  or  (4.6),  one  finds  that  there  are  many 
sets  of  parameters  (X^.,  ,  Z^,  V)  that  result  in  roughly  the  same  D(t). 

Indeed,  the  differences  in  D(t)  between  the  case  of  a  target  moving  slowly 
close  by  a  sensor  and  the  case  of  the  target  moving  more  quickly  further  from 
the  sensor  are  subtle,  and  are  small  compared  to  the  variance  in  estimating 
D(t). 


This  ambiguity  shows  up  In  the  eigen- structure  of  the  matrix  used  in 
finding  the  track  parameter  estimates  via  (4.6).  For  typical  values  of  the 
track  parameters,  the  matrix. 


has  an  eigenvalue  which  is  o{10**-4)  smaller  than  the  rest.  The  coefficeint 
array.  A,  will  therefore  be  very  sensitive  to  any  noise  in  D(t). 


The  eigenvector  corresponding  to  the  troublesome  eigenvalue  is  roughly. 


Therefore,  specification  of  V  and  an  appropiate  modification  of  (4.6)  will 
result  in  track  parameter  estimates  that  are  relatively  insensitive  to  noise 
in  the  measurements  of  D(t).  With  V  specified,  the  modifications  to  (4.6) 
are, 

A  =  minif  -  fj2,  f  =  <&A  (4.9) 


A 


i 


! 


A 


C2D(t1)2t1  C2D(t1)2 

•  • 

'*«*■>**»  a<v2 

(CDItj))4  -  4V2(CD(t1)t1)2 

( CD ( t  ) )4  -  4V2(CD( t  )t  )2 
m  mm 


1 
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The  least  squares  solution  for  A  is  then  given  by. 


*  T  -1  T 
A  =  (<t>  $)  $  f 

and  the  track  parameter  estimates  are. 


XT 

-a./8V 

VT 

3 

(a2/4  +  (a1/8V)2  -  Z2  -  a3/16Z2)1/2 

_ZT  J 

jS  /4ZS 

(4.10) 


(4.11) 
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4.2  TRANSFORM  APPROACH 


In  this  approach,  a  transformation,  which  is  a  function  of  a  depth 
estimate,  is  applied  to  D ( t) .  When  the  depth  estimate  is  accurate,  D(t)  is 
transformed  into  a  second  order  polynomial  whose  coeffeceints  are  easily 
f i tted . 

The  transformation  is  given  by, 

.  4Z.  , 

T  =  (CD  i  ZT  (4.12) 

using  (2.1)  for  D(t), 

T  -  -^|^-[(XT-Vt)2+Y^+(ZT±ZS)2]  +  (XT-Vt)2+Y^(ZTTZ$)2] 

+  i^li[((XT-VT)2+Y^(ZT+Zs)2)((YT-Vt)2+Y^+(ZT-Zs)2j]1/2  (4.13) 

where 

<x  =  zT/zT 

When  a  ~  1  ,  i.e.  Zy  ~  Zj  ,  (4.13)  reduces  to, 

T  =  (XT-Vt)2  +  Y2  +  <ZT  ±  Zs)2  +  O(a-l)  (4.14) 

Fitting  a  polynomial  to  (4.14),  by  least  squares,  for  Instance,  gives, 

V  =  /a^ 

XT  =  -  -i—  a,  (4.15) 

1  2/a^  L 

Yy  =  (a3  -  (ZT  ±  Z$)2  -  X2  )1/2 

where  a^^  ,  e^,  a3  are  chosen  to  minimize 

aaj^t2  +  a2t  +  a3  -  T„  (D( t) ) i 2  , 

ZT 


and  Zj  is  picked  a  priori. 

The  optimal  value  of  Zy  .  and  therefore  (Xy,  Yy,  V)  ,  can  be  found  by 
preforming  a  line  search  over  Zy  ,  using  the  following  as  possible  error 
criteria: 

error  =  UUj t?+a2t.j+a3)  -  TA  (D ( t1 ) ) ]2  (4.16) 

1  ZT 

or,  error  =  jf(D(t)  -  D(t))2, 
i 

where  D( t)  is  given  by  (2.1)  using,  (Xy,  Yy,  Zy,  V)  . 

SENSITIVITY  ISSUES 

As  expected,  the  error,  (4.16)  achieves  a  minimum  only  for  those  cases 
where  D ( t)  is  known  very  precisely  (one  part  in  ten  thousand  for  typical 
values  of  the  track  parameters).  Since  V  can  be  reliably  determined  from 

A 

other  measurements,  a  line  search  over  Zy  can  be  peformed  using, 

error  *  (V(Zy)  -  V)2  (4.17) 

instead  of  (4.16)  as  an  error  criterion.  As  the  transformation  T  results  in 
a  one  to  one  relationship  between  Zy  and  V,  there  are  no  uniqueness 
problems  in  the  determination  of  Zy  •  Knowledge  of  V,  therefore,  allows  use 
of  (4.12)  and  (4.17)  in  estimating  the  remaining  track  parameters  from  noisy 
measurements  of  D ( t) . 

4.3  TWO  SENSOR  EXTENSIONS 

In  this  section,  methods  for  adapting  the  above  techniques  to  inter¬ 
sensor  measurements  in  two  sensor  arrays  are  presented.  In  addition,  methods 
for  combining  Information  from  two  sensors  to  resolve  ambiguities  in 
determining  the  track  parameters  are  developed. 


MODIFICATION  FOR  INTER-SENSOR  DELAY 


It  Is  fairly  straightforward  to  modify  (4.6)  or  (4.12)  so  that  they  can 
be  used  with  Inter-sensor  delay  Information  (instead  of  multipath  delay 
information).  If  the  sensors  are  placed  in  a  vertical  array,  no  modification 
of  either  method  Is  needed,  only  estimated  target  depth  needs  to  be  adjusted 
(see  3.18).  With  the  sensors  placed  in  a  horizontal  array,  changes  to  the 
two  methods  must  be  made. 

The  inter- sensor  delay  in  the  case  of  a  horizontal  array  is  still 
written  as, 

D  =  P-Q  (4.18) 


where,  in  this  case, 

P  =  [(XT-vtcose  -  x$)2  +  (YT-Vtsine)2  +  (ZT-Z$)2] 1/2 

2  1/2 
Q  »  [(XT-Vtcose  +  xsr  +  (YT-Vts1ne)  +  (ZT-Z$)  ]  ' 

and 

(CD)4  +  (cDt)2(-4V2)  +  8(cD)2t  .  (VXT cose  +  VYTsine) 

2  222  22222 
+  (CD)  .  (-Y)(Xy+  X^  +  Yj  +  (ZS-ZT)  )  +  t  .(16X‘V  COS  9) 

+  t  .  (-16X^XtVcos9)  +  (16X^X2)  =  0 

With  the  exception  that  there  are  six  variables  to  be  determined  instead  of 
four,  solving  for  the  coefficeints  of  (4.18)  is  equivalent  to  solving  for  the 
coefficeints  of  (4.2).  The  track  parameters  in  (4.18)  can  be  determined  by 
the  method  described  in  (4.6).  Additionally,  any  a  priori  knowledge  of  track 
parameters  can  be  incorporated  in  a  manner  similiar  to  (4.9). 


The  transformation. 


T  -  (CD  *  hr 


(4.19) 


will  again  yeild  a  polynomial  for  the  correct  choice  of  y  •  The 
coefficeints  of  the  polynomial  can  then  be  fitted,  determining  the  remaining 
track  parameters.  Unfortunately,  the  optimal  choice  of  y  is. 


y  =  2X$(Xt  -  Vtcose) 


(4.20) 


and  the  search  for  y  is  now  over  a  two  dimensional  surface,  which,  without 
a  good  error  criterion,  will  not  yield  accurate  or  computationally  efficient 
results. 

COMBINING  ESTIMATES  FROM  TWO  SENSORS 

Even  though  there  are  many  sets  of  track  parameters  resulting  in 
essentially  the  same  multipath  delay  curve  at  a  sensor,  these  same  sets  of 
track  parameters  will  not  necessarily  generate  similiar  delay  curves  at 
another  sensor  location.  Therefore,  if  information  from  two  sensors  Is 
available,  ambigiuties  in  track  parameter  estimates  can  be  resolved  by 
looking  for  track  parameter  estimates  which  are  consistent  at  the  two  sensor 
locations. 

If  data  is  available  from  two  sensors  placed  in  a  vertical  array,  sets 
of  estimates,  parameterized  by  V  (or  Zy  ),  can  be  constructed  using  either 
(4.9)  or  (4.12).  A  line  search  can  then  be  performed  over  V  (or  Zy)  for  the 
set  of  track  parameters  minimizing  the  error. 


error  =  [Yy(i,l)  -  YT ( 1 ,2 ) ]^ 

where  Yy(*,m)  is  the  Y  estimate  from  sensor  m  based  on  v  =  V 
or,  equivalently. 


(4.21) 


error  =  ^ rCpAQ (1,1)  -  rcpAQ^’2^ 


(4.22) 


where  rrDAn(ji,m)  is  the  mtn  sensor  CPA  range  estimate  to  the  point  (0,0,0) 


based  on  V  =  V  . 

i 


H 

1 


If  the  data  is  from  two  sensors  placed  in  a  horizontal  array,  sets  of 
estimates,  parameterized  by  V  (or  ZT  ),  based  on  individual  sensor  data  can 
be  generated  by  (4.9)  or  (4.12).  Estimates  of  the  bearing  angle, 
parameterized  by  v  (or  Zy  )  can  be  made  based  on  inter- sensor  data  using 
(4.18).  Track  parameter  estimates  can  then  be  chosen  by  finding  the  set  of 


estimates  producing  the  most  consistent  (r^j, 
S6t  ^CPAl*  ^cPAl ’  ^  minimizing  the  error. 


rCPAl ’ 


A 

a) 


,  i.e.,  the 


error  =  [sine  - 


(4.23) 


4.4  COMPUTER  SIMULATION  RESULTS 


In  this  section,  results  of  the  parameteric  fit  methods  of  track 
parameter  estimation  applied  to  simulated  data  are  discussed.  The  equation 
error  method  and  the  transform  method  are  applied  to  delays  simulated  for  one 
and  two  sensor  arrays. 


DATA 


* 

The  data,  d  was  generated  by  calculating  the  true  D(t),  using  (2.1), 
and  adding  white  gausslan  noise.  Since  D(t)  estimates  tend  to  be  less 
certain  near  CPA  (see  appendix  1),  the  additive  noise  Is  scaled  in  proportion 
to  the  delay  value  in  the  case  of  single  sensor  multipath  delay  curves,  and 
in  the  case  of  2-sensor  horizontal  arrays,  the  noise  is  of  constant  variance. 

Figure  4.1  shows  examples  of  the  simulated  data.  In  Figure  4.1,  the 
solid  line  is  the  true  value  of  the  delay;  the  dots  are  the  data  points 
used.  The  signal  to  noise  ratio  in  these  cases  is  about  20. 

EQUATION  ERROR  METHOD 


The  equation  error  method  was  applied  to  data  generated  using, 


Target:  (5t,  1000,  200) 


Sensor:  (0,  0,  300) 


When  no  additive  noise  is  present,  the  track  parameters  are  estimated  to  one 
part  in  104.  With  a  SNR  of  about  20,  however,  the  estimates  are  off  by 
factors  of  100,  due  to  the  sensitivity  problem  discussed  earlier. 

If  the  value  of  V  is  assumed  known,  the  equation  error  method  can  be 
applied  to  the  noisy  data  with  much  improved  results.  Figure  4.2  shows  track 
parameter  estimates  as  functions  of  V  for  the  equation  error  method  applied 
to  a  data  set  with  a  SNR  of  20.  Note  that  at  the  correct  value  of  V,  the 
track  parameter  estimates  are  close  to  the  true  values. 

Sensitivity 

To  illustrate  the  sensitivity  problems,  two  delay  curves,  generated  using 
different  sets  of  track  parameters  have  been  plotted.  Figure  4.3a  and  4.3b 
show  delay  curves  plotted  using  the  estimated  track  parameters  generated  with 
V=5  and  V=10.  Also  plotted  are  data  sets  with  SNR's  of  20.  Both  curves  fit 
the  data  well.  In  fact,  the  differences  between  the  two  curves  themselves  are 
negligible.  Therefore,  in  this  case  additional  knowledge  (knowledge  of  V,  for 
instance)  is  required  to  make  a  unique  choice  of  track  parameter  estimates. 

TRANSFORM  METHOD 


The  transform  method  was  used  to  generate  the  track  parameter  estimates 

A 

plotted  as  functions  of  Zy  in  Figure  4.4.  Figure  4.4a  shows  estimates  made 
from  noiseless  data.  Estimates  made  from  data  with  SNR=20  are  essentially  the 
same.  The  error,  which  Is  the  norm  of  the  difference  between  the  transformed 

A 

delay  curve  and  the  closest  quadratic,  is  plotted  as  a  function  of  Zy  for 
the  two  cases  in  Figure  4.4b. 

As  functions  of  Zy  ,  the  track  parameter  estimates  for  the  two  cases 
are  very  simillar,  and  assume  values  close  to  the  true  track  parameter  values 

A 

for  the  correct  choice  of  Zy  .  The  error  curves,  however,  indicate  that 
practically  any  noise  In  estimating  D(t)  will  lead  to  ambiguities  in 
determining  track  parameter  estimates.  As  V  is  usually  available  through 


DELAY  ESTIMATED/MEASURED 


Figure  4.3a  Estimated  delay  curve  (solid)  using  V=5  m/s  and  the  resulting 
track  parameters  estiated  obtained  with  the  equation  error 
method.  The  dots  are  the  original  data  set. 


DELAY  EST 1 MATED/ MEASURED 


Figure  4.3b  The  same  as  4.3a  using  V=10  m/s  in  the  equation  error  method. 

The  actual  velocity  for  this  example  was  V* 5  m/s.  Note  the 
similarity  between  the  two  estimates. 


t  .  4 


L.fi 
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V, 


other  means,  the  track  parameters  can  be  determined  by  choosing  Zy  as  the  one 
which  gives  the  correct  value  of  V  and  using  Zy  to  determine  the  remaining 
track  parameters. 


TWO  SENSOR  EXTENSIONS 


When  data  from  two  sensors  are  available,  ambiguities  in  choosing  track 
parameter  estimates  can  sometimes  be  resolved  by  looking  for  consistencies 
between  the  two  sets  of  estimates.  In  a  vertical  array,  delay  information 
alone  can  determine  track  parameter  estimates  when. 


VarD 


Yy  $  0(10.|ZS1-ZS2I)  and  <  .1  • 


and  in  a  horizontal  array  when, 


VarD 


Yy  s  0(10.|XS1-XS2I)  and  <  .1 


Both  methods  were  applied  to  two  sensor  data  satisfying  the  above 
constraints.  Computer  simulation  results  are  given  below. 


Vertical  array 


Both  methods  were  used  to  estimate  track  parameters  from  a  two  sensor 
vertical  array.  For  the  simulations,  the  sensors  were  placed  at. 


Sensor  1:  (0,0,300) 


Sensor  2:  (0,0,100) 


with  a  target  at. 


Target:  (5t,  1000,  200) 


and  SNR=2Q. 


The  equation  error  method  was  used  to  generate  track  parameter  estimates 
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as  functions  of  V  at  each  sensor.  The  error  In  (4.21)  is  plotted  in 
Figure  4.5  for  two  different  simulations.  For  this  SNR  and  yt  ,  the  minimum 
of  the  error  function  will  occur  reasonably  near  the  correct  value  of  V. 

Figure  4.6  shows  track  parameter  estimates  as  functions  of  at  each 
sensor.  These  estimates  were  generated  using  the  transform  method.  Figure 
4.7  shows  CPA  to  (0,0,0)  error,  (4.22)  as  a  function  of  .  Again,  for 
these  values  of  SNR  and  <  0(10*  U^-Z^ I )  »  the  minimum  of  the  error 
curve  will  occur  near  the  true  value  of  Z_  . 


Note  that  In  the  above  simul talons  the  Inter-sensor  delay  measurements 
were  not  used.  Using  this  information  further,  estimates  of  Y-  or  CPA  to 
(0,0,0)  as  functions  of  V  or  Z^  could  be  used  to  lower  the  variance  in 
determining  V  or  Z^  . 

Hori zontal  Array 


Track  parameters  were  estimated  based  on  data  from  a  horizontal  array 
with  sensors  at 


Sensor  Locations:  (±150,0,300) 


In  this  case,  the  track  parameters  were  given  by 


=  200 


Target:  V  =  5,  0  =  »/4,  ^  =  0,  Yy  =  650,  ZJ  = 


The  equation  error  method  was  used  to  estimate  track  parameters  at  each  sensor 
as  functions  of  V.  The  bearing  angle  was  estimated  as  a  function  of  V  from 
inter- sensor  data  also  using  the  equation  error  method.  The  error  criterion, 
(4.23)  was  used  to  determine  the  best  estimate  of  V  (see  Figure  4.8).  This 
error  will  reliably  give  good  estimates  of  V,  and  therefore  the  remaining 
track  parameters,  since  the  inter-sensor  angle  estimate  is  an  Increasing 
function  of  V,  whereas  the  individual  sensor  estimate  of  angle  is  a  fairly 
constant  function  of  V. 


parameter  estimates  as  functions  of  k  =  —  for  sensor 


REMARKS 


In  the  case  of  a  single  sensor  in  the  presence  of  a  multipath  reflection, 
both  the  equation  error  method  and  the  transform  method  yeild  accurate 
results  if  the  velocity  or  depth  are  known  Independently.  In  the  case  of  a 
two  sensor  array,  depth  or  velocity  need  not  be  known  a  priori,  but  can  be 
chosen  using  consistency  criteria. 
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5.  CONCLUSION 


In  this  report,  techniques  were  developed  for  extracting  track  parameter 
information  using  sono-bouy  data.  In  the  case  of  a  single  sensor  and  a  target 
moving  in  the  presence  of  a  multipath  reflection,  it  was  shown  that 
differential  delay  and  doppler  information  could  be  obtained  and  used  to  find 
reasonably  accurate  estimates  of  the  depth,  range,  and  velocity  of  the 
target.  In  the  case  of  two  sensors,  additional  range  and  velocity  estimates 
can  be  obtained  from  intersensor  differential  delay  information.  These 
estimates  can  be  combined  to  form  lower  variance  estimates  of  these 
parameters.  Additionally,  depending  on  the  sensor  locations,  further 
estimates  of  depth  and/or  estimates  of  bearing  angle  may  be  found. 


As  a  final  note,  it  is  worth  mentioning  that  a 
spent  developing  software  to  process  sonobouy  data, 
included  correlogram-based  TOOA  estimators  —  SCOT, 
ADEC  linetracking  algorithm.  Listings  are  Included 
data  example  is  presented  in  appendix  2. 


significant  effort  was 
The  software  developed 
PHAT,  and  ML  [4],  and  the 
in  appendix  4.  A  real 
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APPENDIX  1 


DIFFERENTIAL  OOPPLER  EFFECTS  ON  CORRELATION 
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Appendix  1:  DIFFERENTIAL  DOPPLER  EFFECTS  ON  CORRELATION 

If  two  versions  of  the  same  signal  have  relative  doppler  shifts,  the  peak 
of  the  resulting  cross  correlation  will  be  reduced  compared  to  that  of  the 
cross  correlation  of  the  signals  without  a  relative  doppler  shift.  This 
decorrelation  causes  fading  of  correlogram  lines  near  CPA,  where  differential 
doppler  is  greatest. 

To  examine  this  effect,  assume  X(t)  is  a  zero  mean,  gaussian  bandpass 
process,  with  bandwidth  B,  and  center  frequency  Fq  .  The  autocorrel ation  of 
X  is  given  by 


V'1  =  Ro  cos  2" 


(Al.l) 


The  mean  correlation  output,  v  ,  of  the  signal  X(t)  with  itself  is  given  by 

T 

„  7 


w  -  E{^  £  X(t)X(t)dt} 


n  =  R. 


(A1.2) 


The  mean  correlator  output,  M  of  the  signal  X(t)  with  its  doppler  shifted 
version  is  given  by 


T 

2 


U  ■  E{y  £  X(t)X(at)dt) 
which  is  not  in  general  as  large  as  Rg  . 


(A1.3) 


U  _  1  r 

RT‘  .»'l|i-«l  > 
U  0 


ttBT  1 1-al 


uBT| l-a| 


[sinc( 


l)-sinc(-^B-T41~-°1(^0-  1))]  '  (A1.5) 


Plots  of  mean  correlator  output  as  a  function  of  doppler  shift  appear  in 
Fig.  (Al.l).  Figure  (A1.2)  is  a  plot  of  a  typical  doppler  history  for  a 
target  moving  in  a  straight  line  by  an  array.  Several  features  are  easily 
seen  in  Figure  (Al.l).  First,  as  the  relative  doppler  shift  increases,  the 
correlator  output  decreases.  Second,  as  integration  time  increases,  the 
correlator  output  decreases.  Finally,  as  the  bandwidth  of  the  signal 
increases  relative  to  the  center  frequency,  the  correlator  output  increases. 

The  performance  of  the  correlation  method  of  time  delay  estimation  is 
greatly  effected  by  nonzero  differential  doppler.  Fortunately,  due  to  the 
relatively  slow  speeds  of  submarines,  differential  doppler  is  only  likely  to 
cause  problems  except  near  the  time  of  CPA.  See,  for  instance.  Figure  A2.1. 
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Figure  Al.l  Mean  correlator  output  as  a  function  of  doppler  shift,  a 
fg  is  the  center  frequency,  and  g  is  the  bandwidth  of  the 
input  process.  T  is  the  integration  time. 


I 


£ 


„\ 

w  1 


f 


.V 

.  -f 


APPENDIX  2 

REAL  DATA  EXAMPLE 


Appendix  2:  REAL  DATA  EXAMPLE 


This  appendix  presents  an  example  of  track  parameter  estimation  from  a 
real  data  set  taken  from  a  horizontal  array.  The  track  parameter  estimation 
techniques  described  in  Chapter  3  will  be  used.  In  this  case,  e,  V,  1  and 
R  will  be  estimated  from  the  data  and  knowledge  of  the  sensor  locations. 

wr  n 

The  data  was  taken  from  an  array  with  sensors  located  at  (±  136m,  0, 
308m).  Figure  A2.1  shows  R^,  R^»  and  R^  ,  the  correlograms  formed  using 
the  correlation  method  of  time  delay  estimation.  R  and  R  have  been 
energy  normalized  so  that  r^.  (0 )  =  1  V. ;  R  has  been  SCOT  normalized  [4],  As 
the  data  set  was  of  exceptional  quality,  the  differential  delay  information 
needed  for  track  parameter  estimation  can  be  measured  without  further 
processing.  The  immediately  measurable  parameters  are  listed  below: 


D12( | t|  ♦  »)  "  .187  sec 


(A2.1) 


3D1 2 ( t ) 


3t |t:012(t)=0  ^200sec 


,  .5  sec> 
bnncprJ 


=  .0025 


t  -  t  “  35  sec 

rcpa1  rcpa2 


0^(0)  “  .076  sec 


(0 )  "  .074  sec 


-  5.1 


The  relevant  equations  from  the  text  are  as  follows: 


cose  .  CDjltl»H 


(A2.2) 


lT  -  C 

7  p  7Z~ 
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(A2.3) 
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sT  CD(O) 

cdtct  r~ 


( A2.4) 


iA  1 

y  =  3  .  _ i _ 

Cose  tD  -tD 

kcpa1  kcpa2 

Vc  ■  [XsK^Tf  )2-  1]  -  (ZT-Zsl2]'1/2 


(A2.5) 


(A2.6) 


t:019(t)»0 


Using  (A2.2),  one  finds. 


cose  *  1,  e  -  0 


(A2.7) 


Noting  that  the  values  of  0^(0)  and  D2(0)  are  very  close,  and  are  related  to 
RCPA  by  the  ^PA  ran9es  at  sensors  1  and  2  will  be  similar,  confirming 
the  0°  bearing  angle  estimate. 


If  a  0°  bearing  angle  is  assumed,  (A2.2)  gives 


C  -  1450  m/s 


(A2.8) 


with  knowledge  of  C  and  e  ,  V  can  be  estimated  from  (A2.5), 


V  "  7.7  m/s  =  14  KTS 


(A2.9) 


Assuming  ZT  is  constant,  the  two  estimates  of  p  can  be  combined  giving,  from 
(A2.3), 


ZT  -  92m 


(A2.10) 


Using  the  above  value  of  Z^,  R  for  each  of  the  sensors  can  be  calculated 
with  ( A2.4) 


rcpa1  ’  455  m 


(A2.ll) 


R„.  "  470  m 

cpa2 


_V  .a  -V  -i^LSV.  _ V  A  -V.".  AA  .NlVA 


Wil 


The  Y-axis  crossing  can  be  estimated  using  (A2.6), 


Y  -  510  m 
c 


Note  that  since  0 
consistent  with  the  values  of 


°-  Vc  -  RCPA.. 

C 


(A2.12) 


,  and  therefore  the  value  of 

and  R  .  . 

CPA, 


Yc  is 


1  u  2 

From  this  real  data  set,  a  group  of  self-consistent  track  parameter 

values  have  been  obtained  using  various  easily  measured  quantities.  Note  that 

in  this  case,  the  value  of  0  "  0  made  it  difficult  to  determine  Xc  (through 

(55)).  Consequently,  additional  estimates  of  0,  V  and  R„„.  were  not 

CPA 

obtainable. 
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ACCURACY  OF  QEPTH  ESTIMATION  IN  A  MULTIPATH  ENVIRONMENT 
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Palo  Alto,  California  9*304 

ABSTRACT 

The  problem  of  estimating  the  depth  of  an  underwater  source  from  acoustic 
measurements  Is  considered.  The  Cramer-Rao  bound  on  the  variance  of  the  depth 
estimate  Is  evaluated  for  the  cases  of  one  and  two  sensors.  The  effect  of 
multipath  on  estimation  accuracy  is  investigated. 
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1.  INTRODUCTION 


The  problem  of  estimating  the  depth  of  an  underwater  source  Is  of 
considerable  Interest  In  various  surveillance  applications.  Given 
measurements  from  a  vertical  acoustic  array  it  is  possible  to  estimate  the 
angle  of  arrival  of  the  signal.  From  this  estimate  and  knowledge  of  the 
distance  to  the  source  It  Is  possible  to  estimate  its  depth.  If  the  source 
range  Is  not  known  it  can  be  estimated  from  the  curvature  of  the  wavefront  at 
the  array. 

Multipath  propagation  is  a  common  effect  in  the  transmission  of  acoustic 
waves  In  water.  The  presence  of  multipath  propagation  Is  usually  considered 
to  be  an  undesirable  effect  which  complicates  the  processing  of  the  acoustic 
signals  In  localization  and  detection  problems.  However,  in  some 
circumstances  the  presence  of  multipath  can  be  quite  useful. 

As  an  example  consider  the  case  of  depth  estimation  when  only  a  single 
sensor  is  available.  In  the  absence  of  multipath  the  measurements  provided  by 
the  sensor  contain  no  Information  about  the  location  of  the  source.  However, 
if  multipath  Is  present,  the  delay  between  the  signals  propagating  in  the 
direct  and  the  secondary  paths  contain  information  about  the  depth  of  the 
source.  In  the  presence  of  multipath  it  is,  therefore,  possible  to  estimate 
source  depth,  while  in  its  absence  this  is  not  possible. 

In  this  paper  we  study  the  accuracy  of  source  depth  estimation  both  in 
the  presence  and  in  the  absence  of  multipath,  for  the  case  of  one  and  two 
sensors.  These  results  can  be  extended  in  a  straightforward  manner  to  an 
arbitrary  number  of  sensors. 

The  Cramer-Rao  bound  (CRB)  is  used  to  evaluate  estimation  accuracy.  The 
CRB  specifies  the  best  possible  accuracy  achievable  by  any  unbiased  estimator 
Cl].  In  the  case  of  stationary  Gaussian  processes,  the  CRB  can  be  evaluated 
by  a  simple  frequency  domain  formula  as  the  number  of  observations  becomes 
large  [2].  This  simplified  formula  has  been  applied  to  the  estimation  of 
inter-sensor  time  delay,  and  single-sensor  multipath  delay  [3].  The  purpose 
of  this  study  is  to  extend  previous  work  to  include  multipath  effects  in  the 
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two-sensor  case.  The  main  question  is  how  much  better  can  depth  be  estimated 
when  both  multipath  and  inter-sensor  time  delay  Information  are  used  by  the 
estimator. 

Four  cases  are  defined  to  show  a  progression  from  the  single-sensor  case 
to  the  two- sensor  case: 

(1)  Single-sensor  case  with  multipath.  The  multipath  delay  is  used  to 
estimate  depth. 

(2)  Two-sensor  case,  no  multipath.  The  inter-sensor  time  delay  is  used 
to  estimate  depth. 

(3)  Two- sensor  case  with  multipath.  In  this  case  both  inter-sensor  delay 
and  multipath  information  are  used  to  estimate  depth. 

(4)  Two-sensor  case  with  multipath  and  unknown  range.  When  range  is  not 
known,  the  combination  of  inter-sensor  delay  and  multipath  delays 
suffice  to  determine  both  range  and  depth.  The  CRB  for  this  case  is 
also  presented. 

Examples  are  provided  in  each  case  to  illustrate  the  performance  bounds 
as  function  of  system  parameters. 
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2.  ESTIMATING  OEPTH  FROM  SINGLE-SENSOR  MULTIPATH  INFORMATION 


The  singl e-sensor  multipath  case  is  depicted  in  Figure  I.  The  received 
signal  s(t)  can  be  expressed  as 

sit)  *  x(t)  +  gxit-Oj)  +  e(t)  (1) 

where  x(t)  is  the  signal  received  by  direct  propagation  from  the  source  to  the 
sensor,  e(t)  is  uncorrelated  measurement  noise,  and  gx(t-O^)  is  the 
attenuated  and  delayed  version  of  the  source  signal  arising  from  a  surface 
bounce.  The  correspondi ng  power  spectral  density  (PSD)  function  at  the  sensor 
is  given  by 

SsU)  «  SXU)  1  +  ge'jwDl  2  +  Se(cu) 

«  SxU)C(l+g2)  +  2gcosUDx)]  +  SeU),  (2) 

for  tu  in  the  range  [-W.W] ,  where  W  »  2*B  ,  and  B  is  the  bandwidth  in  Hz. 


In  terms  of  the  geometry  (see  figure  I),  the  multipath  delay  is  given  by 

_  1,-  2  ,  .2.1/2  r  2  ,  ,2-1/2,  m 

Oi-^Cy  +  (z+Zj_)  ]  -ly  +  (z-zJ>)  J  }  »  (3) 

where  c  Is  the  speed  of  sound.  Thus,  assuming  the  range  y  and  sensor-depth  Zj 
are  known,  the  source-depth  z  is  determined  by  the  multipath  delay  Dj* 

Using  Whittle's  formula  [2],  the  CRB  for  the  variance  of  the  depth 
estimate  is  given  by 


><«n  2 


7  *  b  ^(iz1)' 


Jj  l  w  C-2gusinu0,]‘ 

*  '  55  i  S2(.) 


2 

^30^ 


a 


»  CRB(D^)(-^--J  . 


where  M  is  the  number  of  Independent  observation  intervals.  (In  continuous 
time,  the  observation  time  is  defined  as  T=N/2B.  In  discrete  time,  when  3 
equals  half  the  sampling  rate,  N  is  the  number  of  sampled  observations.) 
Equation  (4)  can  be  Interpreted  as  a  generalized  form  of  the  CRB,  see  [4], 

2.1  SPECIAL  CASES 

When  the  PSD  of  x(t)  and  e( t)  are  constant  for  |u,|  <_  W  ,  and  zero  for 
lui I  >  W  ,  with  the  signal-to-noise-ratio  (SNR)  defined  as  Sx/Se  ,  we  have 
the  following  approximations  (see  [3] ,[5]): 


( i )  VO,  »  I,  SNR  «  1 


3  I 


CRB(D,  )  •  Jk— .  -L.  ■  * 

TW' g  SNR  NVr  g4 


g  SNR 


(11)  WDX  »  1,  SNR  »  1 
CRB(D, )  »  , 

1  TwV 

which  is  Independent  of  SNR. 


(ill)  ,WD1  «  l 


CRB(D,)  »  IQl2]f|NR+ll  # - 1  .  (7 

1  TW"3  g^SNR^  I.21TD1 

Note  that  the  CRB  depends  on  the  multi  oath  delay  Dj  for  small  delay- 
bandwidth  product  WDj_ ,  but  not  for  WDi  >>  1. 

When  the  range  to  the  source  is  large  compared  to  the  depths  of  the 
source  and  sensor,  the  derivative  |^—  has  a  simple  form: 


«  .  SL 

30^  3Z 


z,  ^  <  <  y 


Note  that  the  variance  of  the  depth  estimate  is  inversly  proportional  to 


rvr.n 
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2 

(2z/y)  ,  where  (22/y)  Is  approximately  the  angle  between  the  source  and  its 

“reflection" ,  as  viewed  from  the  sensor. 

2.2  EXAMPLES 

Figure  2  shows  a  specific  scenario  for  which  the  CRB  was  evaluated. 
Figures  3,  4,  and  5  display  the  normalized  Cramer-Rao  bound,  CR8(z  )/z2  as  a 
function  of  SNR,  for  delay-bandwidth  products  BO^  3  0.1,  1,  and  10, 
respectively.  For  the  geometry  in  Figure  2,  these  time-bandwidth  products 
correspond  to  bandwidths  B  *  3.78,  37.78,  and  377.78,  respectively.  Each 
figure  shows  three  different  values  of  g:  g  3  -0.1,  -0.5,  and  -0.90.  As 
expected,  increased  reflection  magnitude  |g|  Improves  the  ability  to  estimate 
z,  as  does  increased  delay-bandwidth  product. 


3.  ESTIMATING  OEPTH  FROM  INTER-SENSOR  DELAY 


Figure  5  shows  the  general  case  of  two  sensors  configured  in  a  vertical 
array  at  deptns  z^  and  z 2  ,  respectively.  In  this  section  multipath  is 
assumed  absent.  The  inter-sensor  delay,  i.e.,  the  time-difference  of  arrival 
(TDOA)  between  sensors  1  and  2,  is  defined  as 


*  l.r  2  ,  .2,1/2  r  2  ,  ,2,1/2 

°12  a  clty  +  (Z"21)  1  +  (z"z2  1 

The  respective  received  signals  are  given  by  (analogous  to  (U), 


(9) 


Sjft)  -  x(t)  +  e^t)  , 

s2(t)  *  x( t-0 12)  +  e2 ( t)  .  (11) 


where  e,(t)  and  e_(t)  are  assumed  to  be  uncorrelated.  The  power  spectral 

1  ^  *r 

density  of  the  2-vector  process  s(t)  *  [s^ ( t) ,  s2(t)]  is  equal  to 


$(u>) 


j  (uD,  n 

S  (o»)  +  S  (w)  S  (u)e 

X  X 

S  ((u)e’j“°12  SU)  *  Sa  U) 
x  x 


(12) 


For  the  vector  case,  the  asymptotic  form  of  the  CRB  for  the  covariance  of 
the  parameter  estimates  a  becomes  [2], 

Cov(9)  >  J'1  ,  J  *  [Jjj]  , 

JfJ  *  XT  jfw1'*'  ^  1*.  .  <13' 


where  tr ( A )  denotes  the  trace  of  the  matrix  A.  The  CRB  for  0  is  then  given 
by 


CBB(01Z>  •  i?  /u  & 


W  2  2 

T  /  u  SNR  . 

77  /y  U2SNP  <3u 


and  the  CRB  for  depth  z  is  given  by 


CRB(z)  =•  CR8(D12)(!^j)2 


3.1  SPECIAL  CASES 


If  the  noise  and  signal  PSD  functions  are  assumed  flat  in  [-W,W]  and  zero 
elsewhere,  then 

CRB(D12)  »  H-  - -1^-  •  (16^ 

12  ITT  SNR2 

For  SNR  «  1  and  WD  »  1  ,  we  have  (comparing  with  (5)) 

g2  .  ,17) 

CRB(D1) 

At  long  ranges,  y  >>  z,  z^ ,  z^  ,  we  have  (cf.  (8)) 


>012  Zj-lj 


Thus, 


CRB (z)[ two-sensors,  no  multipath] 
wu'mf*..  vajt.rw..'  -w..  rnnvpR-AMI  Ol  Uil  • 


le-sensor,  mu 


(Z2-Zil 


When  g»0  or  z*0,  no  depth  information  is  obtained  in  the  single-sensor 
multipath  case,  while  for  z^  3  ,  no  depth  information  is  available  in  the 

TDOA  case.  These  facts  are,  of  course,  to  be  expected. 
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3.2  EXAMPLES 


Figure  7  shows  the  normalized  CRB,  CRB(z)/z2,  as  a  function  of  SNR  for 
three  different  signal  bandwidths  B  »  3.78,  37.78,  and  377.78.  The  geometry 
is  the  same  as  in  the  single-sensor  case  (figure  2)  with  the  second  sensor 
being  located  at  z^  *  400  .  The  three  curves  may  be  compared  with  the 
results  in  figures  3,  4,  and  5,  respectively.  Figure  8  shows  the  effect  of 
varying  the  geometry,  using  B  3  377.8  Hz  and  z^  3  220,  400,  and  600. 


4.  ESTIMATING  DEPTH  FROM  BOTH  MULTIPATH  AND  INTER-SENSOR  QELAYS 


Fiqure  9  shows  the  geometry  for  the  two-sensor  case  in  which  both 
multipath  and  inter-sensor  delays  are  to  be  estimated.  The  two  multipath 
delays  are  given  by 


.  lrr  2  ,  A  ,2,1/2  r  2 ,,  ,2 ,1/2, 

Di  ■  ^[y  +<z+2i>  ]  -  (y  ]  }  • 


-  l,r  2  .  A  ,2,1/2  ,  2  ,  ,2,1/2, 

°2  +*z+z2>  ]  "  [y  +(z_z2)  J 


and  the  TDOA  is  again  equal  to  (cf.  (10)) 


n  lfr  2  ,  ,2,1/2  r  2  .  ,2,1/2. 

°12  *  ctt*  +lz~h]  1  '  [y  +<z'z2)  1  f 


The  received  signals  at  sensors  1  and  2  are  given  by 


y^(t)  *  x( t-x11 )  +  gx(t-t21)  +  e^t)  , 


y2(t)  »  x(t-x12)  +  gx(t-x22)  +  ez(t)  , 


where 


°1  =  T21  *  TU  ’ 


°2  *  t12  '  t22  ’ 


°12  *  T11  "  t12  ' 


The  PSD  matrix  is  (cf.  (2)  and  (12)) 


Sx(l+gz^2gcos«D1) 


Sxe  ii:(l+ge  i)(l+ge  c) 


S  (id)  *  juiO.«  ju»D.  ju)D5  - 

Sxe  u(l+ge  L)(l+ge  £)  Sx(l+g*>2gcosa)D2) 


The  CRB  for  depth  estimation  is  then 


.  (2 


(28) 


-jwO-,  z-z,  -juQ,  -juD,  z+z,  , 

"( 1+^2a  )(-T7J  -  9i*  k  • 

11  21  c 

4.1  EXAMPLES 


(i)  Different  Bandwidths,  8=3.78,  37.78,  377.78 

Figures  10,  11,  and  12  show  the  normalizede  CRB,  CRB(z)/z2,  for  the  three 
delay-bandwidth  cases  seen  In  the  previous  examples.  The  case  of  g»0  in  each 
figure  is  precisely  the  CRB  curve  obtained  in  the  previous  section  for  TOOA 
based  depth  estimation  (no  multipath).  Mote  that  use  of  multipath  delay 
estimates  can  improve  the  relative  variance  in  the  depth  estimate  by  10  to  15 
dB  when  the  multipath  is  strong. 


(11)  Changing  Geometry 


figures  13,  14,  12,  15  show  a  sequence  of  CRB  curves  In  which  the  second 
sensor  Is  taken  deeper  and  deeper,  I.e.,  z2  =  200,  220,  400,  600  .  All  are 
at  the  large  bandwidth  3  =  377.78.  In  figure  13,  the  g*0  curve  is  at  infinity 
because  for  z ^  *  z2  and  g  *  0  there  is  no  TDOA  or  multipath  information  on 
which  to  base  a  depth  estimate.  We  find  that  increasing  the  depth  of  sensor  2 

Improves  the  bound.  There  Is  more  Improvement  near  g=0  than  near  g*-0.9 
indicating  that  it  Is  the  TDOA  component  that  is  benefiting  from  Increased 
z2  .  This  Is  reasonable  In  view  of  equations  (8)  and  (18)  which  show  that  at 
long  range,  the  derivative  of  z  with  respect  to  TDOA  depends  strongly  on 
z2-z1  while  the  sensitivity  to  multipath  delay  does  not. 
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5.  ESTIMATING  OEPTH  AND  RANGE  FR(M  MULTIPATH  ANO  INTER-SENSOR  DELAYS 


In  the  single-sensor  case.  It  has  been  shown  [6]  that,  assuming  the  range 
error  is  uncorrelated  with  the  multipath  delay  error,  the  variance  of  the 
depth  estimate  becomes  a  sum  of  two  terms  involving  the  respective  variances 
of  range  and  delay: 

Var(z)  -  (|^-]2Var(D1)  +  (|j)2Var(y)  ,  (29) 

where  the  sensitivity  of  y  to  z  is  defined  for  fixed  D^  : 

3Z  a  (30) 

At  long  range  (z,  z^  <<  y)  ,  we  may  use  (8)  to  obtain 

Var(z)  .  Var(0l)  ^  Var(y)  (31) 


For  two  sensors  with  no  multipath  (TDOA  only),  a  similar  situation  is 
obtained. 

In  the  case  of  two  sensors  with  multipath  It  is  possible  to  estimate  both 
range  and  depth.  An  interesting  question  Is  the  following:  how  much  is  depth 
estimation  accuracy  effected  by  the  lack  of  knowledge  of  the  range?  To  answer 
this  question  we  evaluate  the  CRB  on  the  covariance  matrix  of  the  estimate  of 
the  parameter  vector  [z,y]T.  Using  Whittee's  formula  we  note  that 


CRB{[z,y]T}  *  j 


yz  yy 


where  (cf.  (25)), 


■  h  /"Is1”''1  ^  l4*  • 


( 33a) 


(33b) 


Jzy  -Jyz-ir  nsu  , 

1  *  ^  'cl  \“1  ’  f  S  I  \  3S(  m  )  ij 

Jyy  47J.  ls(u,)  1y  jts'u)  • 


(33c) 


The  matrices  S(m)  and  are  defined  in  equations  ( 25 ) - ( 26 ) ,  while 


3  S  ( ui ) 

sy 


3$^  ( at) 
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3S21  (u») 

LT^y— 


3  S  x  2  (  ui ) 


where 


9S^  (01) 

“ay 


952(111) 
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(34) 


(35) 


-JmD,  -JuD.  -  -JmD.  , 

+  g2e  Z(1+g,e  l)-*-- (i*g2«  Z)  — - 

z  1  t22  i  T11 


-JmD,  -JmD2  , 

-  V  1 


Exampl es 


Figures  16,  17,  and  18  show  the  normalized  CRB,  CRB (z)/z^,  for  the  three 
bandwidths  considered  in  the  previous  examples  ( BD^ *0.1,  1,  10)  .  In  this 
particular  geometry,  the  TDOA  component  of  the  deoth  estimate  is  independent 
of  the  range  due  to  the  hyperbola  of  constant  TDOA  degenerating  to  a  straight 
line.  Consequently,  for  g»0,  the  CRB  curves  reduce  to  the  TDOA-only  case  of 
figure  7.  Note  that  at  small  30  ^  ,  the  presence  of  multipath  only  makes  the 
estimates  worse,  while  at  intermediate  and  high  BO^  ,  the  depth  estimate  is 
improved  by  multipath. 


fi 

i 

Figure  19  shows  the  results  at  high  BD^  for  a  different  geometry 
uj|  (2^*600  instead  of  400).  In  this  case  the  TDOA  comoonent  of  the  depth 

estimate  depends  on  range.  This  figure  compares  with  figure  18.  In  the  non- 
degenerate  geometry,  the  unknown  range  adds  significantly  to  the  uncertainty 
i n  depth . 


s 

I 


5.  OISCUSSIOM 


The  performance  analysis  presented  in  this  paper  shows  the  potential 
usefulness  of  multipath  information  in  enhancing  the  accuracy  of  depth 
estimation.  The  results  presented  here  are  asymptotic,  and  need  to  be 
verified  by  finite  data  simulations.  Two  related  issues  of  practical  interest 
are: 

(I)  The  Case  of  Unknown  g 

The  earlier  analysis  assumed  perfect  knowledge  of  the  relfection 
coefficient  g.  When  g  is  unknown,  it  can  be  estimated  from  the 
data.  In  [6]  we  have  shown  for  the  examples  considered  here  that  the 
asymptotic  accuracy  of  depth  estimation  is  not  affected  by  the 
need  to  estimate  g,  provided  that  ’4)  »  1,  where 
D  *  maxtO^  02>  D^}. 

(II)  The  Case  of  Unknown  SNR 

Similarly  to  the  case  of  unknown  g,  it  can  be  shown  [61  that  the 
need  to  estimate  SNR  does  not  affect  the  asymptotic  depth  estimation 
accuracy,  provided  that  WD  »  1. 

Finally,  we  note  that  these  results  can  be  extended  in  a  straightforward 
manner  to  the  case  of  M  sensors.  In  this  case  SU)  and  $(u)  will  be  Mrti 
matrices.  The  entries  of  these  matrices  can  be  easily  evaluated  by  obvious 
extrapolation  from  the  case  of  M«2. 
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Figure  1:  Single  sensor  geometry 
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Figure  2:  Single  Sensor  Example 
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Figure  3:  Single  sensor  with  multipath,  B  =  3.78 
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Figure  7:  Two  sensors  no  multipath  example  —  varying  bandwidth 
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Figure  11:  Two  sensors  with  multipath  --  23  3  400.  B  =  37.77 
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Figure  16:  Unknown  range  —  Z2  *  400,  B  *  3.78 
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Figure  19:  Unknown  range  —  z 2  *  600,  B  *  377.70 
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APPENDIX  -  SOFTWARE  LISTIN6S 

The  software  developed  to  compute  and  display  the  CRB  for 
the  various  test  cases  was  written  in  the  CtrlC  language. 

CtrlC  is  a  software  product  of  Systems  Control  Technology. 

The  CtrlC  code  below  is  contained  on  the  SCT  VAX  directory 
NETVAX:  :  dra2:  Cmts.  bound,  cr.  depth 3*.  *.  While  connected  to  that 


B 

£ 

U 

a 

i 


directory/  the  monitor  command  &CMCM  will  run  an  example  test 
case.  The  command  SCMCM  invokes  the  VAX/VMS  DCL  command  file 
CNCH.COM  which  is  also  listed  below.  A  summary  of  the  various 


a 


files  is  given  first. 

cmcm.  com  -  Run  CtrlC  and  feed  it  emcm.  ctr 

cmcm.  ctr  -  Outer  loop  of  CtrlC  commands  for  different  SNR/  BTz  etc 

tdoa.  ctr  -  Compute  CRB  for  case  of  two-channel  TDOA 

scmp.  ctr  -  Compute  CRB  for  case  of  single-channel  multipath  data 

mcmp.  ctr  -  Compute  CRB  for  case  of  multichannel  multipath  data 

mctd. ctr  -  Compute  multipath  delays  and  tdoa  from  geometry 

mpsn.  ctr  -  Adjust  SNR  to  eliminate  signal  boost  due  to  multipath 

simp,  ctr  -  Simpson's  rule  for  numerical  integration 

The  file  cmcm.com  loads  the  CtrlC  function  library 
dra2:  CMTS.  ML.  SPLIB1SP.  LIB  which  contains  many  simple  utility  functions 
not  implemented  directly  in  CtrlC.  The  contents  of  SP. LIB  are 
described  in  CSmithS4b3. 


b 
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CMCM.  COM  -  Initiation  Of  Test  Cases 


*!  cmcm.com  -  Co  into  CtrlC/  load  SPLIB/  and  DO  cmcm.  CTR  to  obtain  the 
<!  crb's  for  the  1-channel  and  2-channel  multipath  situation 

*! 

*  run  dral: Csctlib.  ctrlclctrlc.  exe/nodebug 
d iary  >cmcm.  log; 

lib  DRA2:  CMTS-  ML.  SPLIBlspi 
j  *  sqrt(-l); 

//  echo»l; 

//  term* 'vt  12 
term*'ptx 
hard  *  'ptx 
redh  >cmcm.  ptx 
do  cmcm; 
exit; 

*  pra2  cmcm. ptx 
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*  pra2  ctrlc.  ptx 

*  pra2  cmcm.  log 
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CMCM.  CTR  -  Outer  Loop  Of  CtrlC  Commands  For  Different  SNR,  8T, 
Etc. 


//  C 3*cmcm< ) 
deff  mcmp 
deff  scmp 
deff  tdoa 
deff  mpsn 
deff  mctd 
deff  simp 
dotdoa»0; 
doscmpaO; 
dompsnr*0; 
z  1=*200; 

snr=C-20,  -10,  0,  10,  203;  nsnr*5; 
g«C-0.  001,  -0.  01,  -0.  051,  ng*3; 
vzpn=»ones  <nsnr,  ng+1 ); 

ntd=9;  //  max  number  of  tdoa-only  plots 

vztd=ones (nsnr,  ntd); 

itd-1; 

FOR  z2aC6003... 

FOR  z=*C3003,  .  . 

FOR  y=C30003,  .  . 

FOR  B=»C103/.  026476,  .  . 

Cmpdl,  mpd2,  tdl23«mctd  <  zl,  z2,  z,  y  >»  .  . 
b  t l=mpd 1*B; b t2«mod2*B;  b  t3»ab s ( td 12)*B;  .  . 

strS-C';  BT1» ',  cvf  s  <  btl,  2) ,  '  BT2-',  cvfs(bt2,  2),  '  BT12»',  cvfs(bt3,  2)  3 

Nfrq»max  <  C2*round  <  b ) ,  1003 )  +  l;  .  . 

strO«*C  ';  Nfrqa ',  c vs  < round  <Nfrq)  )  3;  .  . 

disp(Cstr3,  str03); . . 

tmpl«C'  z  1» '»  cvs  (  z  1 )  3;  .  . 

tmp2»C  '  z»',  cv*(z),  '  y«',cvs<y),'  b*'»  cvf  s  <  b,  2)  3;  .  . 

IF  dompsnr*l, tmp2*Ctmp2,  'C '3;  ELSE  and;., 
str  1-C tmp  1,  '  z2*',cvs(z2),  tmp23;  .  . 
str4»Ctmpl,  tmp23;  .  . 
gstr* '  . 

IF  DoTdoa-l,  .  . 

FOR  i«l : nsnr,  .  . 

vzpn  (  i,  ng+1 >«tdoa (  snr  <  i  ),  zl,  z2,  z,  y,  B ); .  . 

IF  itd<*ntd»  vztd < i,  itd >»vzpn< i, ng  +  1 );  ELSE  end;., 
d isp ( C 'TDOA  Var  <  z )■ cvs ( .  5*db ( vzpn( i, ng+1 ) ) ) ,  '  dB'3);.. 
end, . .  i 
ELSE  end; . . 
itd«itd+l; . . 

FOR  ig*l: ng, . . 

tmp«cvf s ( g (  ig  ) ,  1 ) ; . . 
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gstr=Cgstr,  tmp3;  .  . 

IF  ig<ng,  gstr*Cgstr,  ',  '3; ELSE  end; . . 
sstr-'  .  . 

FOR  i*l :  nsTvr,  .  . 

tsnr*mp sn  <  snr ( i  ) ,  z  1,  z2 ,  z,  y»  3,  g  < ig  )  ) ;  .  . 

IF  domp snr*l, csnr*tsnr; ELSE  c»nr»»nr (  i  ) ;  end;  .  . 
tmp*C  cvs  (  snr  (  i  >  >  <  '  (  cvf  s  <  tsnr,  2)  >  ' )  ' 3;  .  . 
str3*C  '  snr*  '»  tmp  3;  .  . 
sstr*Csstr,  t mp3; . . 

IF  i<nsnr,  sstr*C s str,  ',  ' 3; ELSE  end;.. 
disp(Cstrl,  str2,  str33); . . 

IF  doscmp*l . . . 

vzsc(i<  i g ) *scmp ( c  snr,  zl,z,y,B.g(ig).  Nf  rq  >  i  .  . 
dispCC'lch  Var < z )* ' , cvs (0.  5*db ( vzsc C i, i g ) ) ) ,  '  dB'3);.. 
ELSE  end.  .  . 

C  tmp.  vy  pn,  vy  zpn3*mcmp  ( c  snr.  zl.z2.  z.y.B.  g(ig).  Nfrq) ;  .  . 
vzpn< i.  ig )*tmp;  .  . 

disp<C'2ch  Var  (  z  )*'»  cvs  < 0.  5*db  ( vzpn  <  i.  ig  ))).  '  dB'3);.. 
disp(C'2ch  Var  <  y  )»',  cvs  <0.  5*db  ( vypn  ))»  '  dB'3>;.. 
disp(C'2ch  Xv<zy)*',cvs<0.  5*d  b (abs(vyzp))).  '  dB  '  3  ) . . . 
•nd. .  .  i 
end. .  .  ig 
IF  doscmp*l.  .  . 

•rasa; . . 

p  lot  ( snr, db(vzsc)/2)>.  . 
title(C'lch  Var<  z  ) ',  str43  >;  .  . 
xlab<C'SNR(dB>;  £*'»  gstr, strO. str53>;  .  . 
y lab (  '10 log  <Var <  z ) /z**2> . 

IF  norm<t«rm<  1:  3)-'ptx  ')O0.  replot;  ELSE  and,-.. 

ELSE  end; . . 
erase; . . 

p  lot  < snr,  db  < vzpn) /2) ; . . 
title(C'2ch  Var  <  z )'.  str  13 );.  . 
x  lab  <  C  'SNR <  dB ) ;  G* gstr,  0 ',  strO,  str 5 3  ) ;  .  . 
gstr*  '  . 

y lab ( 'lOlog (Var <  z )/z**2)  '); . . 

IF  norm(term<  1:  3)-'ptx  '  )O0»  replot;  ELSE  and;., 
and, . . B 
end, . . y 
and, . .  z 
and.  //  z2 
IF  DoTdoa*l,  .  . 
erase;  .  . 

p lot<  snr, db(vztd)/2);  .  . 

title <C 'All  TDOA  Var(z);  LAST  ',str43>;.. 
xlab<C  'SNR(dB);  LAST'.  strO,  str53);  .  . 
y lab (  '  lOlog  < Var <  z  > /z**2> . 

IF  norm(  term(  i;  3)-'ptx  '  )O0»  raplot;  ELSE  end;.. 

ELSE  end; 
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TDOA.  CTR  -  Compute  CRB  For  Case  Of  Two-channel  TDOA 


//  Cvzpnl  »  tdoa <  snr.  zl.  z2.  z.  y.  B ) 

// 

//  Compute  CRB  for  source  depth  given  2-channel  TDOA  measurements 

// 

//  zl  *  Sensor  1  depth 
//  z2  ■  Sensor  2  depth 
//  z  »  Target  depth 

//  y  «  Target  range  projected  to  surface  (meters) 

//  B  *  Target  bandwidth  (Hz) 

//  Snr  *  Signal  to  noise  ratio  in  dB  at  each  sensor 

c  =  1500;  //  Speed  of  sound  (m/sec) 

snrl  ■  I0**(snr/10);  snr2  «  10** < snr/ 10 ) ; 

til  *  sqrt(y*y+<  z-zl  )**2)/c; 

tl2  ■  sqrt(y*y+( z-z2)**2)/c; 

d  -  1 12  -  til; 

dztmp  »  z I*tl2-z2*tl l-z*d; 

if  abs ( d ztmp Keps. 

disp('  tdoa:  numerical  failure,  tdoa  not  a  function  of  depth ')» 

vzpn*l.  0;  ..  create  line  at  OdB  .. 

return; 

dzdd  s*  c*c*t  1  l*t  12/d  ztmp; 

vard  ■  (3/<2*pi*B)**2)*(l+snrl*snr2)/(snrl*snr2); 
vzpn  ■  (vard*dzdd**2)/(z*z ); 
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SCMP. CTR  -  Compute  CRB  For  Case  Of  Single— channel  Multipath  Data 

//  Cvzpnl  ■  scmp  (  snr.  z  1.  z.  y .  B»  g.  Nfrq ) 

// 

//  Compute  per-sample  CRB  for  1-channel  multipath  situation 
//  Sampling  rate  is  assumed  to  be  2B  rad/sec 

// 

//  zl  »  Sensor  1  depth 
//  z  ■  Target  depth 

//  y  ■  Target  rang*  projected  to  surface  (meters) 

//  B  *  Target  bandwidth  (Hz) 

//  g  *  Multipath  attenuation 

//  Snr  *  Signal  to  noise  ratio  in  dB 

//  Nfrq  »  Number  of  frequency  samples  uniform  on  CO. BJ 

// 

//  Derived  constants 


c  ■  1500;  cs  ■  c*c; 
f  »  b*CO: Nfrq-1 1/Nfrq; 
w  *  2*pi*f; 

Sx»l;  S*  «  10**<-snr/10) ; 


//  Speed  of  sound  (m/s*c) 

//  Hertz  frequency  axis 
//  Radian  frequency  axis 
//  PSD  in  band  C-B. B3 
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til  *  sqrt < y*y  +  <  z-z  1  )**2) /c; 
t21  *  sqrt  <  y*y  +  (  z  +  z  1  )**2) /c; 

D  *  t2l-tll; 

dDd z  «  ( z+zl )/(cs*t21 )  -  (  z-z 1 )  /  ( c s*tl 1 ) ; 

Sy  *  ( ( lrg**2)*0NES(l.  Nfrq)+2*g*cos(w*D> >*Sx  +  Se*ones ( 1 .  Nf r q ) ; 

SydD  ■  -2*Sx*g*w.  *sin(w*D); 

si  *  SydD.  / <SyrEPS*ONES< 1. Nfrq) )» //  Specific  information  for  delay  tst 
JD  *  SUM< si.  *si ) /<2*Nfrq) ;  //  par-sample  CRB  Var<D>  .  ge.  1/JD 

vzpn  *  ( l/<JD*<dDdz**2)+EPS) )/<z*z ); 


imat— 


MCMP.  V2  -  CRB  For  Multichannel  Multipath  Data.  Known  Range 


CvzpnH  »  mcmp  (  snr,  z  1  ,  z2,  z>  y .  B,  g,  Nfrq  ) 

Compute  CRB  for  2-channel  multipath  situation 


Sensor  1  depth 
Sensor  2  depth 
Target  depth 

Target  range  projected  to  surface  (meters) 
Target  bandwidth  (Hz).  Must  be  from  0  to  .  5 
Common  multipath  attenuation 
*  Signal  to  noise  ratio  in  dB 
*  Number  of  frequency  samples  uniform  on  CO.  B3 


//  y  =»  rarget  range 

//  B  »  Target  bandwi 

//  q  -  Common  multip 

//  Snr  »  Signal  to  n 

//  Nfrq  »  Number  of 
// 

//  Derived  constants 


j  *  sqrt(-l)j 

gl  *  g;  //  Relative 

g2  *  gi  //  Relative 

c  =*  1500;  //  Speed  of 

f  »  b*CO: Nfrq-ll/Nfrq; 
w  ■  2*pi*f; 

Sx  1*1;  Sx2»l; 

snrl  »  snr;  snr2  *  snr; 

Sel  »  10*-*(-snrl/10); 

Se2  ■  10**<-snr2/10); 

//  disp <  'disabling  channel 

//  Multipath  time  delays 


attenuation  of  2ndary  path  to  channel  1 
attenuation  of  2ndary  path  to  channel  2 
sound  (m/sec) 

//  Hertz  frequency  axis 
//  Radian  frequency  axis 
//  PSD  of  signal  in  band  C-B.B3 
//  SNR  in  channels  1  and  2.  resp. 
//  PSD  of  noise  in  band  C-pi.piJ 
//  PSD  of  noise  in  band  C-pi.pil 

1');  Sxl-O; 


•  sqrt ( y*yr< z-zl )**2) /c; 

•  sqrt < y*y+( 1 )**2>/c; 

»  sqrt(y*y+(z-z2)**2)/c; 

■  sqrt <y*y+< z+z2)**2) /c; 

max <C(t22-t 12 >*B.  < t21-t 1 1 ) *81 ) ; 

■  Nfrq/BT;  //  Number  of  integration  points  per  multipath  spectral 


if  ppc<20.  disp(C'  Warning. 
disp(C'  which  means  only 


BT«'»  cvf  s(BT.  2).  '  while  Nfrq» '.  cvs  (Nfrq )  1 
'.  cvfs(ppc.  2).  '  integration  points  per  eye 


cycle 

>;  .  .  :v 

1  e  '  3  ) ;  *£ 


//  derivatives  of  delay  mrt  z 


3 

9 


\\ 


* 

«ij< 

*.■1 


s. 

V 


'v 


'1'.. 


Vi 


CS  *  C*Ci 

tllr  ■  < z-z 1 ) / < cs*tl 1 ) ; 
t21 z  -  < z+zl >/<cs*t21>; 
t I2z  »  < z-z2)/(csetl2>; 
t22z  -  < zez2)/(eset22)i 

//  Derivatives  of  sourc e— to— sensor  spectra  utrt  delays 
//  si/2  «  sl/tll  tll/2  +  sl/t21  t2i/2 

sill  »  -(2*Sx l*g 1 )*(w.  *sin< <tll-t2i)*w) );  //  sl/tll  »  -  sl/t21 

s 1 2  a  si  1 1*( tl 1 2— 121 Z  ) i  //  slll*tll2  +  3 121*t21 2 

//  s2/z  »■  s2/t  12  1 12/ z  *  s2/t22  t22/z 

s212  ■  —  (2*Sx2*g2)  *<ui.  *s in <ui*(  1 12— 122 )  )  ) ;  //  »  -s222 

s2z  -  s212*< t 12z  -  t22z ) * 

//  sl2z  ■  conj(s2l2>  *  sum(  i<  j=*l»  2)  sl2/tij  tij/z 

eju»  »  exp<jeui)> 

ejll  »  exp  <  j*t«*tll  ); 

e j 12  *  exp ( jew»tl2); 

ej21  *  exp ( jew#t21 ) ; 

ej22  »  exp  <  j*u»*t22); 

k I»ejl2+g2*ej22; 

k2»conj<ejll>gl*ej21 ); 

t* jewesqrt (Sx l#Sx2) ; 

xll  »  -t.  econj(ejll).  ekli  //  x  expands  to  “S12/t“ 

x21  »  -g  l*t.  *con  j  <  e  j21 ) .  *k  1; 
x  12  «  t.  *e  j  12.  *k2; 
x22  »  g2*t.  *e  j22.  *k2; 

sl2z  »  x  1  let  1 1  z  *■  x  12*tl2z  *■  x21et21z  *  x22*t22ii 

si  »  Sx l*abs < ones < 1,  Nfrq )eg lee j 1 1.  econj < ej21 )) *e2  +  Sel*ones < 1/ Nfrq) ; 
s2  *  Sx2*abs(ones<l.  Nfrq)eg2*ejl2.  *conj<ej22> )**2  +  Se2*ones ( 1. Nfrq) ; 
sl2  *  sqrt < Sx leSx2)*k2.  ekl; 

//  Do  sisz  »  S** ( — 1 )  S/2 

//  dl  means  d  In  ...  dlij  means  SA-1  dS/dz  Ci*j3 
detS  »  si.  *s2  -  sl2.  *con  j  <  sl2) ; 
dill  -  < s2.  *sl z-sl2.  *con j ( sl2z > ) . /detS; 
d  1 12  »  <s2.  *sl2z-sl2.  *s22).  /detSi 
dl21  *  (  si.  »conj  ( s!2z  )-conj  ( sl2).  *slz ). /detSi 
dl22  *  { si.  *s2z-conj ( sl2>.  *sl2z ).  /detSi 

//  Trace  dl**2 

td  Is  ■  dill,  edlll  e  dl22.  edl22  e  2ed  1 12.  ed  121; 
stdl  *  max < Ceps.  abs( sum< tdls ) ) 3 ) i 

if  stdl»eps.  d isp ( 'numerical  failure,  vanishing  info  matrix'); 
si  ■  sum(tdls( 1:  2:  Nfrq) ); 
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im 


> 

> 


I 


s2  »  sum(tdls(2:  2:  Nfrq) ); 

errv  *  C ( s l-s2) /std 1.  ab * < si +s2 >-s td 1 3; 

if  sum(abs (irrv) )>0.  01. 

d i sp (  '  Crelat ive  half-sum  difference.  03 
errv.  .  . 
else  end; 

tinf  *  std 1/ < 2*Nfrq) ; 
v:pn  *  ( 1/tinf ) / ( z*z ) ; 


& 

4 


(numerical  integration  check):  ') 

Bf 


MCMP.  CTR  -  CRB  For  Multichannel  Multipath  Data.  Unknown  Range 


//  Cvzpn.  vypn.  vyzpl  ■  mcmp  (  snr.  z  1.  x2.  z>  y.  B.  g>  Nfrq ) 

// 

//  Compute  CRB  for  range  and  depth  estimation  given  2-channel 
//  multipath  and  tdoa  measurements. 

// 

//  zl  a  Sensor  1  depth 
//  z2  *  Sensor  2  depth 
//  z  “  Target  depth 

//  y  a  Target  range  projected  to  surface  (meters) 

//  B  =  Target  bandwidth  (Hz).  Must  be  from  0  to  .  5 
//  g  *  Common  multipath  attenuation 
//  Snr  *  Signal  to  noise  ratio  in  dB 

//  Nfrq  a  Number  of  frequency  samples  uniform  on  CO.  83 

// 

//  vzpn  «  Cramer-Rao  lower  bound  for  the  estimation  of  z 

//  vypn  a  Cramer-Rao  lower  bound  for  the  estimation  of  y 

//  vyzp  «  Cramer-Rao  lower  bound  for  cross-variance  of  y  and  z 

// 

//  Derived  constants 


i 

i 


i 


j  a  sqrt(-l); 

gl  *  gj  //  Relative 

g2  *  g;  //  Relative 

c  »  1500;  //  Speed  of 

f  »  b*CO: Nfrq-13/Nfrq; 
w  a  2*pi*f; 

Sz  1*1;  Sx2-1; 

snrl  »  snr;  snr2  »  snr; 

Sel  »  10** (-snr 1/10) ; 

Se2  «  10**(-snr2/10) ; 

//  disp ( 'disabling  channel 


attenuation  of 
attenuation  of 
sound  (m/sec) 

// 

// 

// 

// 

// 

1');  Sx  1*0; 


2ndary  path  to  channel  1  N 
2ndary  path  to  channel  2 

//  Hertz  frequency  axis 
Radian  frequency  axis 
PSD  of  signal  in  band  C-B.B3 
SNR  in  channels  1  and  2.  resp. 
PSD  of  noise  in  band  C-pi.pil 
PSD  of  noise  in  band  C-pi»pi3 


• ') 


2.3 

5 


>5 

••j 


//  Multipath  time  delays 

til  *  sqrt ( y*y*( z-z 1 ) **2) /c; 
t21  »  sqrt ( y*y*< z+z 1 )**2) /c; 
tl2  ■  sqrt ( y*y*( z-z2)**2) /c; 
t22  ■  sqrt ( y*y*< z*z2>**2> /c; 
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BT  *  max < C < t22-tl2)*B.  < t21-t 1 1 > *B 3 > ; 

poc  *  Nfrq/BT;  //  Number  of  integration  points  per  multipath  spectral  cyci 
if  ppc<20>  disp<C'  Warning.  BT* '. cvf s (BT. 2) .  '  while  Nfrq*  '. cvs <Nfrq> 3 ) ; 
disp(C'  which  means  only  '. cvf s ( ppc.  2) .  '  integration  points  per  cycle'3) 

//  derivatives  of  delay  wrt  z 

cs  3  c*c; 

tllz  »  ( z-t 1 > / ( c s*tl 1 > ; 
t21z  -  ( z  +  z 1 > / <  c  s*t21 ) » 
t 12z  -  ( z-z2>/(cs*tl2>i 
t22z  *  <z«-z2)/<  cs*t22); 

//  derivatives  of  delay  wrt  y 

c  s  *  c  *c  ; 

1 1 1 y  *  y / ( c  s*t 1 1 )  i 
t21y  *  y/(cs*t21); 
tl2y  *  y/(cs*tl2>. 
t22y  *  y/(cs*t22); 

//  Derivatives  of  source-to-sensor  spectra  wrt  delay,  depth,  and  range 


//  sl/z  »  sl/tll  tll/z  +  sl/t21  t21/ z 
sill  =»  -<2*Sx  l*g  1 )  ■*(  w.  *s  in  <  <  1 1  l-t21  >  *w)  ) ; 
slz  *  si 1 1*( tl 1 z-t21 z ) ; 
sly  »  slll*<tlly-t2ly ); 

//  s 2/ z  -  s2/t 12  t 12/ z  +  s2/t22  t22/z 
s212  »  -<2*Sx2*g2)*(w.  *sin<w*(tl2-t22> > >; 
s2z  a  s212*< t 12z  -  t22z ) < 
s2y  a  s212*(tl2y  -  t22y ); 


//  Sl/tll  a  -  Sl/t21 
//  s 1 1 1 *t 1 1 z  +  Sl21*t21z 
//  s 1 1 l*tl ly  +  sl21*t2ly 

//  a  -s222 


//  sl2z  «  conj(s21z)  *  sum< 

ejw  a  e x p ( j * w > ; 

ejll  ■  exp < j*w*tll ); 

ejl2  *  exp ( j*w*tl2) ; 

ej21  ■  exp < j*w*t21 ); 

ej22  »  exp < j*w*t22); 

k l*e  j 12rg2*e  j22> 

k2»c  on  j <  e j 1 1+g l*e j21 ) j 

t»  j  *w*s  qr  t  <  S  x 1 *S  x  2 ) ; 

Ill  »  -t.  *c on j < e j 1 1 ) .  *k 1 j 
x21  *  -gl*t.  *conj(ej21).  *k  1 
x  12  ■  t.  *e jl2.  *k2; 
x22  a  g 2*t .  *ej22.  *k2; 
sl2z  »  x  1  l*t  1 1  z  x  12*t  12z 

sl2y  a  x l l»t 1 ly  r  xl2*tl2y 


sum<i.j«1.2)  sl2/tij  tij/z 


expands  to  "S12/t’‘ 


x21*t21z 

x2l*t21y 


x22*t22z; 
x22*t22y ; 


si  »  Sx l*ab s ( ones ( 1.  Nfrq)+g l*e j 1 1.  *con j < e j21 ) ) **2  +  Se i*ones < 1 . Nf r q ) ; 
s2  «  Sx2*abs(ones(  1.  Nfrq  >+g2*e  j  12.  *con  j<ej22)  )**2  *■  Se2*ones  ( 1.  Nfrq )  i 
sl2  a  sqrt(Sxl»'}x2)*k2.  *kl; 


iCHin 


detS  »  si.  **2  -  sl2.  *con  j  <  sl2)  ; 

//  Do  sisz  *  S**(-l)  S/z.  Below#  dzij  means  S^-l  dS/dz  Ci#j3 
dzll  ■  (  s2.  *sl  z-sl2.  *con  j  <  sl2z  )  ) .  /detS; 
d z  12  *  (s2.  *sl2z-sl2.  *s2z).  /detS; 
d:21  *  ( si.  *con j ( sl2z )-con j ( sl2) .  *sl z ) .  /detS. 
dz22  =  ( si.  *s2z-con j ( sl2) .  *s 12z ) .  /detS; 

//  Trace  dz**2 

tdzs  »  dzll.  *dzll  +  dz22.  *dz22  ♦  2*dzl2.  *dz21; 

stdz  =»  simp(tdzs)# 

ctdz  *»  simp  (  td  z  s  ( 1 :  2:  Nfrq  )  ) ; 

disp('  Jll:  Csimp ( n >-simp ( < n-1 ) /2) 3/ simp ( n ) 3  =  '); 

(stdz-ctdz)/max(Cstdz#  ep s3 ) 

disp('  Jll:  C  simp  ( n  ) -s  imp  (<  n-1 ) /2  >  3  *  '); 

std  z-c  td  z 

if  abs ( imag < std z ) )>eps#  d i sp < 'real i ty  failure/  depth  part');  stc 
stdz  *  real <  stdz  ); 

Jmtx  *  0NES(2>;  //  Allocate  Fisher  information  matrix 

Jmtx(l#l)  =  stdz/2; 

//  Do  sisy  *  S**(-l)  S/y.  Below#  dyij  means  S'*-l  dS/dy  Ci#j3 
dyll  *  (  s2.  *sly-sl2.  *con  j  (  sl2y  )  ).  /detS# 
dy  12  *  (  s2.  *sl2y-sl2.  *s2y  ).  /detS# 
dy21  »  (si.  *con j ( s I2y )— c on j ( sl2 ) .  *sly).  /detS# 
dy22  *  (  si.  *s2y-con  j  (  sl2) .  *sl2y  ) .  /detS; 

//  Trace  dy*#2 

tdys  *  dyli.edyll  *  dy22.  *dy22  *  2*d y  12.  *d y21 ; 

stdy  ■  simp(tdys); 

ctdy  »  simp ( tdys (1: 2: Nfrq) ); 

disp('  J22:  Csimp <n )-simp << n-1 ) /2) 3/simp ( n ) 3  *  '); 

< stdy -ctdy) /max ( Cstdy#  ep s3 ) 

disp('  J22:  Csimp (n )-simp ( (n-1 ) /2> 3  *  '); 

stdy-ctdy 

if  abs( imag ( stdy ) )>eps#  disp ( 'reality  failure#  range  part');  std 
stdy  *  real < stdy ) ; 

Jmtx (2# 2)  »  stdy/2; 

//  Trace  dz*dy 

tzys  «  dzll.edyll  +  dzl2.  *dy21  +  dz21.  *dyl2  +  dz22.  *dy22; 

stzy  ■  simp (tzys); 

ctzy  »  simp ( tzys< 1: 2: Nfrq) ); 

disp('  J12:  Csimp (n)-simp ( (n-1 )/2) 3/simp (n) 3  *  '); 

(stzy -ctzy) /max (Cabs(stzy)# ep  s  3 ) 

disp('  J12:  Csimp (n)-simp < (n-1 ) /2) 3  ■  '); 

stzy-ctzy 

if  abs < imag (stzy ) )>eps#  disp ( 'real ity  failure#  range  part'); 
stzy  *  real (stzy); 

if  abs ( stzy Xeps.  disp('mcmp:  zero  off-diagonal  in  info  matrix') 
Jmtx(1.2)  *  stzy/2# 

Jmtx  (2,  1)  *  Jmtx(l#2); 
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long;  joit* 
short; 

d i sp ( ' c ond i t ion  number  of  J  matrix: ');  CQND(Jmtx) 


Jinv  *  INV(Jmtx); 
v  z  p  n  »  Jinv(l.l)/(z*z); 
vypn  »  Ji nv ( 2. 2) / ( y*y ) ; 
vyzp  =  jinv ( 1.  2) / ( y*z ) ; 


MCTD.  CTR  -  Compute  Multipath  Delays  And  Tdoa  From  Geometry 


a  inctd(zl>  z2i  z>  y) 

compute  multipath  delays  in  each  channel 

to  surface  (meters) 


//  Cmpd 1,  mpd2. tdl23 

// 

//  Given  geometry. 

//  and  tdoa  between  channels. 

// 

//  zl  *  Sensor  1  depth 

//  z  a  Target  depth 

//  y  *  Target  range  projected 

// 

C  =  1500;  CS  a  c*c; 

til  a  sqrt ( y *y  +  ( z-z 1 )**2) /c; 

t2l  a  sqrt < y*y+< z+:l )**2> /c; 

mpDl  a  t21-tl 1; 

tl2  a  sqrt < y*yr( z-z2)**2) / c; 

t22  *  sqrt ( y*y+< z+z2)**2) /c; 

mpD2  »  t22-t 12; 

td 12  a  tl2-t 1 1; 


//  Speed  of  sound  (m/sec) 

//  source  to  sensor  1  direct 
//  source  to  sensor  1  reflected 
//  multipath  delay,  sensor  1 
//  source  to  sensor  2  direct 
//  source  to  sensor  2  reflected 
//  multipath  delay,  sensor  2 
//  tdoa.  channel  1  minus  2 


MPSN.  CTR  -  Adjust  SNR  To  Eliminate  Signal  Boost  Due  To  Multipath 


//  Ccsnrl  a  mp sn ( snr .  z 1 .  z2.  z .  y . B.  g  ) 

// 

//  Compute  SNR  correction  due  to  signal  power  gain 
//  in  the  presence  of  multipath. 

// 

//  zi  a  Sensor  1  depth 
//  z  *  Target  depth 

//  y  ■  Target  range  projected  to  surface  (meters) 

//  B  ■  Target  bandwidth  (Hz) 

//  g  a  Multipath  attenuation 

//  Snr  a  Signal  to  noise  ratio  in  dB  given  no  multipath 

// 

//  Derived  constants 


c  *  1300;  cs  =  c*c; 

W  »  2*p  i*B; 

til  »  sqrt  (  y  *y  +  <  z-z  1  >**2 ) /c ; 
t21  =  sqrt ( y *y+< z+z 1 )**2) /c; 

01  a  t21 -t 1 1 ; 

UIDl  »  w*Dl; 

facl  a  i+g**2+2*g*SIN<WDl ) /WD1; 

tl2  *  sqrt (y*y+( z-z2)**2) /c; 
t22  =  sqrt  <  y*y  +  <  z  +  z2)  *-*2) /c; 

02  »  t22-tl2* 

WD2  a  Ui*02i 

fac2  =  i+g**2+2*g*SIN(  WD2> /WD2; 

f  snr*-lO.  if  abs(g)a0. l*  disp< 
CB*D1* B*D2, B*ab  s ( tl2-tll ) 1 

fac  a  sqrt < f ac l*f ac2) * 
if  facCeps*  fac*eps;  disp('mpsn: 

csnr  »  snr/fac; 


//  Speed  of  sound  <m/sec) 

//  Radian  band  limit 

//  source  to  sensor  1  direct 
//  source  to  sensor  1  reflected 
//  multipath  delay*  sensor  1 


//  source  to  sensor  2  direct 
//  source  to  sensor  2  reflected 
//  multipath  delay*  sensor  2 


'BT  for  MP1.MP2*  and  TDOA:  '), 


//  I  don't  know  if  this 
numerical  failure:  zero  delay'); 

//  corrected  SNR 


is  be 


SIMP.  CTR  -  Simpson's  Rule  For  Numerical  Integration 


//  C  inti  “  simp ( f ) 

//  integrate  f  using  Simpson's  rule.  Integral  is  NORMALIZED, 
n  ■  ma  *  <  s  i  z  e  <  f  )  ) ; 

if  round (n/2)*2an.  disp('simp:  Need  odd  number  of  points'); 
if  n<3>  disp('simp:  Need  at  least  5  data  points'); 
int  »  ( f ( 1 )+f (n )+4#sum< f (2:  2:  n-1 ) >+2*sum< f (3:  2:  n-2) ) ) /3; 
int  »  int/ ( n-1 ) » 
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APPENDIX  4 


CORRELOGRAM  GENERATION  PROGRAM 


The  following  are  listings  of  the  correlogram  and  line  tracking  software 


-yJr. 


,  -* ..’•"v** -k„"' _ 

.  «v^Viiv.rMi<v^-  AJV-V.  AVJif-' 


DRA2 : [ ABEL ] XRSDA . FOR ; 7 


24-FEB- 1986  12:21 


S 


i 


I 


c 


c 


i 


I 


s 

’•  »- 


c 


C 


uc 
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options/check-all 
program  xrsda 

XRSDA  cakes  . i2  data  from  xfile  and  yfile  and  produces  auto- 
and  cross-  correlograms  and  spectrograms  stored  as  .  r4  in 
rxyfile,  sxyfile,  rxxfile,  etc. 

The  ith  line  of  the  cross -correlogram  is  +/-1024  lags  of  the 
biased/unbiased  ML/PHAT/SCOT/unormalized  cross -correlation  estimate 
based  on  nspl*nrps*256  points,  beginning  at  the  (d*i)th  point 
of  x  and  y.  The  lines  of  the  auto-correlograms  are  1024  lags 
of  the  auto-correlation  estimate. 

Correlograms  are  computed  as  if ft{ spectrograms } .  The  spectrograms 
are  computed  as  peroidograms  based  on  nspl  non- overlapping  averages 
of  nrps*256  point  8k  fft's  of  x  and  y. 


XRSDA  'decimates'  the  correlograms  and  spectrograms  for  display  on 
the  DeAnza.  r42mas,  daform  and  fix  should  be  used  to  create  plot 
files  on  floppy  disks . 


declarations 


implicit  none  ! !  input/computation  variables 

character*20  xfile , yfile , rxyfile , sxyfile 

character*20  rxxfile , sxxf ile , ryyf ile , syyf ile 

character*16  nrmlzn 

integer  i , ii , iii , j ,k, 1 , itemp 

integer  nspl ,nrps , sdf ,m,d,nl,nrpl, frr 

integer  lpco,hpco,nrm 

integer*2  x(256 , 2048) ,y(256 , 2048) 

real*4  r (2048 ,3) ,s(512,3) , rtemp(3) , scl(3) 

logical  gtyes ,xey ,bce 


real*4  temp (1024*64)  ! !  **** 

integer  a,b,c  !!  **** 


integer  status, ierr  !!  ap  variables 

integer  xadr , xxadr , yadr , yyadr , wadr , sadr 

integer  xyadr , tempadr , f ts , nc ,ne 

real*4  apbufr (8192 , 3) , apbufs(4097 , 3) , apbufw(8192) 


integer  xch.ych, sxych, rxych 


! !  i/o  variables 
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integer  sxxch , rxxch , syych , ryych 


initialization 


xfile— 'bvdata:mpcl. i2 *  !!  set  default  input  values 

yfile-'bvdata:mpc2 . i2' 

rxyfile-' eg: rxy . r4' 

sxyfile-' eg: sxy . r4' 

rxxfile— ’ eg: rxx. r4 ’ 

sxxfile-' eg: sxx. r4' 

ryyf ile- ' eg : ryy . r4 ' 

syyf ile- '  eg :  syy .  r4 ' 

xey-. false . 

nrmlzn-' nonePHATSCOTML  ' 

nl-512 
nrps-16 
m— 512 
sdf-8 
d-1 

nspl-4 
lpco— 4096 
hpco-1 
frr-1 
nrm-0 

xadr-1000  ! !  ap  defaults 

xxadr-9200 
yadr-17400 
yyadr-25600 
xyadr-33800 
tempadr-42000 
sadr-50200 
wadr-55000 

nc— 4097 
fts-8192 

a-0 

b-1024*64-l 
c-2 

xch-1  ! !  i/o  defaults 

ych-2 
rxych-3 
sxych-4 
rxxch-13 
sxxch-14 
ryych-23 
syych-24 


input 


!  !  **** 
!  i  **** 
!  i  **** 


S 


0 


call  clrscrn 
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call  otstr( ' Program  XRSDA 
call  otstr(' 
call  otstr(' 


computes  the  atuo-  and  cross-' ,1) 
correlogram  and  spectrogram' , 1) 
of  x,y.  Output  for  DeAnza',5) 


call  gtstr(' Enter  the  name  of  the  input  file  for  x  ',xfile,20> 
call  gtstr('Enter  the  name  of  the  input  file  for  y  ’,yfile,20) 
call  gtstr('Enter  the  output  file-name  for  Rxy ' , rxyfile , 20) 
call  gtstr('Enter  the  output  file-name  for  Sxy ' , sxyfile , 20) 
call  gtstr(' Enter  the  output  file  name  for  Rxx ' , rxxf ile , 20) 
call  gtstr(' Enter  the  output  file  name  for  Sxx' , sxxf ile , 20) 
call  gtstr('Enter  the  output  file  name  for  Ryy ' , ryyfile , 20) 
call  gtstr('Enter  the  output  file  name  for  Syy ' , syyfile , 20) 
call  otstr('  ' ,2) 

call  gtint( 'Number  of  records  per  processing  segment ', nrps , 16 , 2) 
itemp-int(2048 ./float(nrps) ) 

call  gtint( 'Number  of  segments  per  correlogram  line ' ,nspl , itemp , 1) 
call  gtint( 'Number  of  records  skipped  between  lines d, 99999 , 1) 
call  gtint( ' Number  of  the  first  record  read’ , frr, 99999 , 1) 
call  gtint( ' Number  of  lines  computed' , nl , 512 , 1) 
call  otstr ( '  ’ , 2) 

call  gtint('0-no  normalization,  1-PHAT,  2-SCOT,  3-ML' ,nrm, 3 , 0) 
bce-gtyes( 'Biased  correlation  estimates  (y/n) : ' ) 
call  otstr( '  ' , 2) 

call  gtint('Enter  low  pass  cut  off  --  4096  bins  ' , lpco ,4097 , 1) 
call  gtlntC Enter  high  pass  cut  off  --  4096  bins ' ,hpco , lpco , 0) 


%  0 

,f.;50 

'£41 

^2 

143 

$4, 

iSso 

170 

">• 

.180 
"  * 

V94 

■ 


call  clrscrn  ! I  print  input  values 

call  otstr('Data  entered: ’,2) 

write(6,120)  xfile 

format('  data  source  l:',a20) 

write(6,130)  yfile 

format ('  data  source  2:',a20) 

write(6,140)  rxyfile 

format('  Rxy:',a20) 

write(6,150)  sxyfile 

formate  Sxy:  '  ,a20) 

write(6,141)  rxxfile 

formate  Rxx:',a20) 

write (6, 142)  sxxf ile 

formate  Sxx:',a20) 

write(6,143)  ryyfile 

formate  Ryy:',a20) 

write(6,144)  syyfile 

format('  Syy:',a20) 

call  otstr ( ' S ' ,1) 

write(6,160)  (nspl*nrps) 

formate  #  of  256  pt  records  per  correlogram  line  :',i6) 
write(6,170)  d 

format('  #  of  256  pt  records  skipped  between  lines :',i6) 
write(6 , 196)  (frr-1) 

formate  #  of  records  skipped  before  the  1st  line  :',i6) 
write(6,180)  nl 

format('  #  of  correllograra  lines  : ' , i4) 

write(6,194)  nrps 

formate  #  of  records  processed  per  8k  fft:',i4) 
call  otstr( ' ,1) 

write (6 , 192)  (nrmlzn(nrm*4+l : nrm*4+4) ) 
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92  formate  power  specturm  normalization:  '  ,a4) 
write(6,193)  bee 

93  formate  biased  correllation  estimates 12) 
call  otstr('b' , 1) 

write(6,195)  hpco/2 , lpco/2 

95  formate  correlogram  passband  --  ( '  , i5 , '  , '  , i5 , ' )  Hz') 
call  otstr( ' c* ,2) 


if  (gtyes(' Values  OK  (y/n) : ' )  .eq.  .false.)  go  to  200 

! !  reenter  if  mistakes 


post- input  initalization 


if  (xfile  .eq.  yfile)  xey-.true.  !!  auto/cross  flag 

ne-nrps*256 

nrpl-nrps*nspl 


do  i-1,8192 

apbufw(i)-0. 

enddo 

do  i-1 , ne - 1 
apbufw(i)-1.0 
apbufw(8193-i)-1.0 
if  (.not. bee)  then 

apbufw ( i ) -float ( ne ) /float ( ne - i+1 ) 
apbufw(8193-i)-float(ne)/float(ne- i) 

endif 

enddo 


open  files 


call  opn (rxych.rxy file, 'new' )  !!  output  files 

call  opn(sxych,sxyfile, 'new' ) 

call  opn(rxxch,rxxfile, 'new' ) 

call  opn(sxxch,sxxfile, 'new' ) 

call  opn ( ryych, ryy file , 'new' ) 

call  opn(syych,syyfile, 'new' ) 


call  opn(xch,xfile, 'old' )  M  open  xfile 

if  (.not.  xey)  call  opn(ych,yf ile , ’ old’ )  !!  open  y--if  necessary 


[ ABEL ]XRSDA. FOR; 7 

main  loop  --  get  data,  compute  correlogram  lines 


do  230  j-1 ,nl 
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if  (j  .eq.  1)  then 


!!  first  data  set  --  read 
! !  nrpl  records 


call  rrecxy(xch,x(l,l) ,ych,y(l,l) , frr , frr+nrpl- 1 , xey) 

!!  not  first  data  set  -- 
if  (d  .ge,  nrpl)  then  !!  d>#  --  fresh  data 

ii—  < j -l)*d+frr 

iii- ii+nrpl-1 

call  rrecxy(xch,x(l , 1) ,ych,y(l , 1) , ii , iii ,xey) 

else  ! !  make  room  for  d 

! !  records  new  data 
do  260  i-l,nrpl-d 
do  270  ii-1,256 
y(ii,i)-y(ii, i+d) 
x(ii,i)-x(ii,i+d) 
continue 

! !  get  d  new  records 
!!  put  at  end  of  x,y 
iii— ( j -l)*d+nrpl+frr 
ii-iii-d+1 
i-nrpl-d+1 

call  rrecxy(xch,x(l, i) ,ych,y(l, i) , ii, iii, xey) 

endif 


endif 


ap  computations 


if  (j.eq.l)  then 

call  apinit(0 ,0 , status) 
call  vclr(0, 1,65535) 
call  apwr 

call  apput(apbufw,wadr,8192,2) 
call  apwd 


endif 


call  vclr(1000, 1,53500) 


! !  initalize  ap 
! !  clear  ap 
!  !  load  constants 


!  !  clear  all  but  window 


do  l-l,nspl 

call  vclr(xadr, l,2*ne) 
call  vclr(yadr , 1, 2*ne) 


!  !  process  x,y 
! !  clear  data  area 
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call  apwr  !!  x,y  -->  ap 

call  apput(x(l , (l-l)*nrps+l) ,xadr,ne,l) 
call  appuc(y(l, (l-l)*nrps+l) ,yadr ,ne  ,  1) 
call  apwd 


call  vflt(xadr , 1 ,xadr , 1 ,ne) 
call  vflt(yadr , 1 ,yadr , l,ne) 


! !  convert  to  reals 


call  rfft(xadr,fts,l) 
call  rfft(yadr,fts,l) 
call  rfftsc(xadr , fts , 3 , 1) 
call  rfftsc(yadr,fts,3,l) 


! !  take  8k  fft' s 
! !  scale/undo  wierd  packing 


call  vmov(xyadr , 1 , tempadr , 1 , fts+2)  !!  conj (X)*Y+sum  -->  xyadr 

call  cvma(xadr , 2 ,yadr , 2 , tempadr , 2 , xyadr , 2 ,nc, -1) 


call  vmov(xxadr , 1 , tempadr , 1 ,nc)  !!  conj (X)*X+sum  -->  xxadr 

call  scjma(xadr , 2 , tempadr, 1 , xxadr , 1 ,nc) 


call  vmov(yyadr, 1, tempadr, l,nc)  !!  conj (Y)*Y+sum  -->  yyadr 

call  scjma(yadr , 2 , tempadr , 1 .yyadr , 1 ,nc) 


enddo 


! !  get  new  x,y 


if  (nrm.eq.2)  then 


! !  SCOT 


call  vmul (yyadr , 1 .xxadr , 1 , tempadr , 1 ,nc) 
call  vsqrt(tempadr , 1 , sadr , 1 ,nc)  !!  (xx*yy)**0 . 5 

call  crvdiv(xyadr , 2 , sadr , 1, tempadr , 2 ,nc)  !!  SCOT 

call  vmov( tempadr , 1 , xyadr , 1 , 2*nc) 


endif 


if  (nrm.eq.l)  then 


! !  PHAT 


call  vmov(xyadr, 1 , tempadr , 1 , nc*2) 

call  cvmags( tempadr , 2 .xyadr , 1 , nc)  !!  magnitude**2 
call  vsqrt (xyadr , 1 , sadr , 1 , nc)  !!  sqrt(mag**2) 
call  crvdiv( tempadr , 2 , sadr , 1 , xyadr , 2 , nc)  !!  PHAT 


endif 


if  (nrm.eq.3)  then 


!  !  ML 


V'  DRA2 : 

if 
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call  cvmags(xyadr,2,sadr,l,nc) 
call  vsqrt(sadr , 1 , tempadr , 1 ,nc) 
call  vclr(xadr,l,2*nc) 

call  crvdiv(xyadr , 2 , tempadr , 1 , xadr , 2 , nc ) 
call  vclr(yadr,l,2*nc) 
call  vmul (xxadr , 1 , yyadr , 1 , yadr , 1 , nc ) 
call  vsub(sadr , 1 , yadr , 1 , xyadr , 1 ,nc) 
call  vdiv(xyadr , 1 , sadr , 1 , tempadr , 1 ,nc) 
call  crvmul(xadr , 2 , tempadr , 1 , xyadr , 2 , nc) 
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! !  mag**2(xy) 

! !  mag(xy) 

! !  phase (xy) 

! !  xx*vy 
! !  xx*yy-xy**2 
! !  xy**2/( . . .) 
!  !  ML 


endif 


call  vmov(xyadr , 1 , tempadr, 1 , f ts+2) 
call  cvmags ( tempadr , 2 , sadr , 1 ,nc) 
call  vmov(sadr, 1 , tempadr , 1 ,nc) 
call  vsqrt( tempadr , 1 , sadr , 1 ,nc) 


! !  mag  of  Sxy 


call  apwr 

call  apget(apbufs(l , 1) , sadr ,nc ,2) 
call  apget(apbufs(l , 2) , xxadr ,nc ,2) 
call  apget ( apbuf s (1,3) ,yyadr,nc,2) 
call  apwd 


xy- ->apbuf 
xx- ->apbuf 
yy- ->apbuf 


if  (hpco.ne.O)  then 

call  vclr (xyadr, l,2*hpco) 
call  vclr (xxadr, l.hpco) 
call  vclr (yyadr , 1 ,hpco) 

endif 


if  (lpco.ne.4097)  then 

call  vclr(xyadr+2*lpco, 1, 2*nc-2*lpco) 
call  vclr(xxadr+lpco,l,nc-lpco) 
call  vclr(yyadr+lpco , 1 ,nc-lpco) 

endif 


! !  Ip  filter 
! !  Sxy,  Sxx,  Syy 


! !  hp  filter 
!!  Sxy,  Sxx,  Syy 


call  vmov(xxadr , 1 , tempadr , 1 , nc) 

call  vclr (xxadr , 1 , 2*nc) 

call  vmov( tempadr , 1 , xxadr , 2 ,nc) 


! !  make  Sxx  cplx 


call  vmov (yyadr , 1 , tempadr , 1 ,nc) 

call  vclr(yyadr,l,2*nc) 

call  vmov( tempadr , 1 .yyadr , 2 , nc) 


! !make  Syy  cplx 


call  rf ftsc(xyadr , fts , -3 , 0) 
call  rfftsc(xxadr,fts, -3,0) 


! !  pack  for  ifft 


WW 
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c 

c 

c 


c 


c 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


320 

c 

c 

c 

c 


c 


call  rfftsc(yyadr,fts, -3,0) 


call  rfftb(xyadr , tempadr , fts , -1) 
call  vmul ( tempadr , 1 , wadr , 1 , xyadr , 1 , f ts ) 


!  !  (un)biased  ifft 
!  !  in  xyadr 


call  rfftb(xxadr, tempadr , fts , -1) 
call  vmul( tempadr , 1 , wadr , 1 , xxadr , 1 , fts) 


!  !  (un)biased  ifft 
!  !  in  xxadr 


call  rfftb(yyadr, tempadr , fts , -1) 
call  vmul ( tempadr , 1 , wadr , 1 , yyadr , 1 , f ts ) 


! !  (un)biased  ifft 
! !  in  yyadr 


call  apwr 

call  apget(apbufr(l , 1) , xyadr , fts , 2) 
call  apget(apbufr(l,2) .xxadr , fts , 2) 
call  apget(apbufr(l , 3) , yyadr , fts , 2) 
call  apwd 


!  xyadr- ->apbufr 
!  xxadr- ->apbufr 
!  yyadr- ->apbufr 


if  (j.eq.l)  then 
do  k-1 , 3 
rtemp(k)-0. 
enddo 

scl(l)-sqrt(apbufr(l,2)*apbufr(l, 3) ) 
if  (nrm)  scl(l)-float(fts) 
scl(2)-apbufr(l,2) 
scl(3)-apbufr(l , 3) 

endif 


!  !  initalize  output 


1-0  ! !  decimate  S 

do  320  i-ltnc-l 
do  k-1, 3 

rtemp (k)-rtemp (k) +apbuf s ( i , k) 
enddo 

if  (mod(i,8) .eq.O)  then 
1-1+1 
do  k-1, 3 

3 ( 1 , k)— 10*logl0( 1 . e- 30+rtemp(k)*f loat ( f ts) / ( 8 . *scl (k) ) ) 

rtemp (k)-0 . 

enddo 

endif 

continue 


1-0 

do  k-1, 3 
rtemp (k)-0 . 
enddo 


136 


v. 


bJ 


'A 


i 


JW 

a 

js 

Vi 

ts 


s 

•c, 

ii 


5k 

i* 


**  >  A  ,N  . 


i  a: 


DRA2 : [ABEL]XRSDA. FOR; 7 
do  321  i-1,512 


!  !  fill  Rxx.Ryy 


24-FEB-1986  12:21 


do  k-2,3 

r(i,k)-apbufr(i ,k)/scl(k) 

rtemp(k)-0 . 

enddo 

21  continue 


1-0 

do  k-1,3 
rtemp(k)-0 . 
enddo 

do  322  i-1,1024  !!  fill  Rxy 

rtemp(l)— rtemp(l)+apbufr ( i , 1)**2 
rtemp(2)— rtemp(2)+apbufr (7168+i , 1)**2 

if  (mod(i,4) .eq.O)  then 
1-1+1 

r ( 256+1 , l)-rtemp(l)/( (scl(l)**2)*4. ) 
r(l , 1 ) — r temp (2)/((scl(l) **2 ) *4 . ) 
rtemp(l)-0 . 
rtemp(2)-0. 

endif 

12  continue 


call  wrec(rxych, r(l , 1) , 4*j - 3 ,4*j )  !!  write  r,s 

call  wrec ( sxych , s ( 1 , 1) , 4*j -  3 , 4*j ) 

call  wrec(rxxch, r(l , 2) ,4*j -3 ,4*j ) 

call  wrec(sxxch, s(l , 2) ,4*j -3,4*j ) 

call  wrec(ryych, r(l , 3) ,4*j -3 ,4*j ) 

call  wrec(syych, s(l , 3) ,4*j - 3 ,4*j ) 


call  otstr ( ' : ' , 0) 


0  continue 


call  cls(xch)  • !  close  data  files 

call  cls(ych) 


call  cls(rxych) 


!  !  close  output  files 


DRA2 : [ ABEL ] XRSDA . FOR ; 7 


24-FEB-1986  12:21 


call  cls(sxych) 
call  cls(rxxch) 
call  cls(sxxch) 
call  cls(ryych) 
call  cls(syych) 

call  otstr ( '  ' , 1) 

call  otstr( 1  ...  files  written', 1) 


stop 

end 


include  ' dra2 : fabel . jsalibjj saio . for/list  *  !!  get  i/o 
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file  i/o  subroutines  --  records 
written  as  512  byte  blocks 


subroutine  opn(ch,nam, oon) 

opens  file  w/  name-nam  on  unit  ch 

implicit  none 
integer  ch.ios 
character*(*)  nam 
character*(3)  oon 


! !  open  file 


open(unit-ch,name— nam, status-oon, err-100 , iostat-ios , 

recordtype-' fixed' , recordsize-512/4 , access-' direct' , 
organization-' sequential' ,blocksize-512*64) 
return 

call  prterr('opn  error’ ,ios)  !!  print  error 

return 

end 


subroutine  cls(ch) 

closes  file  on  unit  ch 

implicit  none 
integer  ch.ios 

close (unit-ch , status— ' keep ' , err-100) 
return 

call  prterr('cls  error', ios) 

return 

end 


! !  close  file 


! !  print  error 


subroutine  rrecxy(chx, xbuf , chy ,ybuf , sr , fr ,xey) 


reads  starting  sr  finnishing  fr:  chx-->xbuf,  chy-->ybuf 
if  xey  sets  ybuf-xbuf  --no  chy  exsists 


implicit  none 

integer  chx, chy , nrcds , sr , fr , i 

byte  xbuf (abs((fr-sr+l) *512) ) , ybuf (abs( (fr-sr+l)*512) ) 
logical  xey 
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c 

c 

call  rrec(chx , xbuf , sr , fr)  !!  read  x's  records 

c 

if  (xey)  then 

nrcds-fr-sr+1 
do  200  i-1 , nrcds*512 
200  ybuf (i)-xbuf (i) 

else 

call  rrec(chy,ybuf , sr , fr) 

endif 
c 

return 
end 
c 
c 
c 
c 
c 

subroutine  rrec(ch,buf , sr , fr) 
c 

c  reads  fr-sr  512  byte  records  from  ch  to  buf,  starting  at 

c  record  sr  finishing  at  fr 

c 

implicit  none 

integer  sr , fr , ch, ios , i , s , j 
byte  buf (abs( (fr-sr+l)*512) ) , temp (512) 
c 
c 

if  (sr  .gt.  fr)  then  !!  starting  record  after 

ios-999  ! !  ending  record 

goto  200 

endif 
c 

do  100  i-sr.fr  !!  read  records  sr...fr 

c 

read(unit-ch , rec— i , err-200 , iostat-ios ) temp 
c 

s-(i-sr)*512  !!  put  read  buffer  in 

do  j -1,512  !!  return  variable,  buf 

buf (j+s)-temp(j )  !!  512  bytes/record 

enddo 
c 

100  continue 

c 

return 
c 

200  call  prterr('rrec  err', ios)  !!  print  errors 

return 
end 
c 
c 
c 
c 
c 

subroutine  wrec(ch,buf ,sr, fr) 
c 

c  writes  buf  to  ch  in  512  byte  records  -- 


!!  xey-->set  y-x 


!!  else  read  y's  records 
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■"‘c  starting  at  record  sr,  finnishing  at  fr 


i 


implicit  none 

integer  s , ch , sr , fr , ios , i , j 

byte  buf (abs( (fr-sr+l)*512) ) , temp (5 12) 


9 


* 


c 


c 


^200 

V) 

M 


C 


c 

c 


if  (sr  .gt.  fr)  then 
ios-999 
go  to  200 

endif 

do  100  i-sr,fr 


!  !  starting  record  after 
!  !  last  record 


!!  write  records  sr...fr 


s-(i-sr)*512 
do  j-1,512 

temp( j )-buf ( j+s)  !!  fill  write  buffer  -- 

enddo  ! !  512  bytes/record 

write (unit-ch, rec-i , err-200 , ios tat-ios) temp 
return 


call  prterr('wrec  error', ios)  !!  print  errors 

return 

end 


i 


terminal  i/o  routines 


c 


,v 


•fbo 

200 


s 


subroutine  otstr(str ,nlf) 

prints  str,  on  the  current  line,  followed  by  nlf  cr/linefeeds 


implicit  none 
integer  nlf,l 
character*(*)  str 

write(6 , 100)str  !!  prints  string  to  terminal 

format(lh+,a,$) 

if  (nlf  .It.  1)  return  !!  prints  lf's 

do  200  1-1, nlf 

write(6 , 300) 

format( '  ' ) 

continue 

return 

end 
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c 

c 

c 

c 


c 


100 

c 

200 

c 


300 

c 

c 

c 

c 

c 

c 

c 

c 


c 


100 

c 

200 

400 

c 

c 

c 


subroutine  gtstr(pt,var,n) 

gets  a  string  from  the  terminal  --  prompts  w/  pt,  var=string 

var  has  length  n<101 


implicit  none 
character*(*)  pt,var 
integer  n,i, nchrs 
character*100  invar 

i-min(len(var) , index(var , '  ')) 

write(6,100)  pt.var  !!  write  prompt,  default 

format('  ',a,'  (’ ,a<i>, ’ ) : ' ,$) 

read(5 , 200 )nchrs ,( invar ( i : i) , i-1 , nchrs)  !!  read  srting 
format(q , lOOal) 

if  (nchrs  .eq.  0)  return  ! !  return  default 

do  300  i-l,n 
var ( i : i ) - '  ' 

if  (i  . le .  nchrs)  var (i : i)-invar ( i : i)  !!  fill  return  variable 

return 

end 


subroutine  gtint(pt, igr,mx,mn) 

gets  integer  from  terminal  --  prompts  user  w/  pt,  accepts 

igr  in  range  (mn,mx) 


implicit  none 
character*10  itstr 

integer  igr , mx , mn , it , nchrs , i , lmn , lmx 
character*(*)pt 


i-2+max(l , int(logl0(abs ( float (igr) )+ . 1) ) ) 

write(6 , 100)pt , igr  ! !  write  prompt,  default 

formate  ',a,'  ('  ,  i<i>, ' ) :  '  ,  $) 

read(5 ,400)nchrs , (itstr(i: i) , i-1, 10)  !!  get  string 

format(q , 10a!) 


if  (nchrs  .eq.  0)  return 
call  s2i( it , itstr) 


!  take  default  if  no 
!  entry 

!  convert  str-->int 


if  ((it  .le.  mx)  .and.  (it  . ge .  mn))  then  !!  check  range 

igr-it 
return 

endif 

lmn-2+max(l, int(loglO(abs(floatj (mn) )+. 1) ) ) 
lmx-2+max(l , int(loglO(abs(floatj (mx) )+. 1) ) ) 
write (6 , 300)mn,mx 


!  !  not  in  range 
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formate  The  (min, max) are 

go  to  200 

end 


c 


( ' , i<lmn> , ' , ' , i<lmx> , ' ) 
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reenter: ' , $) 


'-’t 


c 


»  * 

~c 


c 


.  n 

,v 


$ 

ft 


C-Qoo 

-Jbo 


c 


subroutine  s2i(itr,str) 

converts  string  --  str  --to 

implicit  none 
integer  a,l 

integer  itr , i , j , temp ,pt 
character*(*)  str 

temp— 0 
pt-1 

l-len(str) 

integer  --  itr. 

do  100  i-0,1-1 
a— ichar (str(l-i:l-i) ) 

! !  itr-sum{ str(i)*10**i) 

if  (a  .eq.  48)  go  to  200 

!  a-'0’ 

if  (a  .eq.  45)  go  to  300 

!  a-'-' 

if  ((a  .It.  49)  .or.  (a  .gt. 

temp-temp+pt* ( a- 48 ) 

pt-pt*10 

continue 

itr- temp 

return 

itr-- temp 
return 

end 

57))  go  to  100 

!  a  not  in 

subroutine  gtrel(pt,rel,mx,mn) 

gets  real*4  from  terminal  --  prompts  user  w/  pt,  accepts 

rel  in  range  (mn,mx) 


implicit  none 
character*10  rlstr 
real*4  rel,mx,mn,rl 
integer  nchrs , i , lmx , lmn 
character* (*)pt 

write(6 , 100)pt, rel  !!  write  prompt,  default 

formate  ’.a,’  (', f<5+max(l , int(loglO(abs (rel)+. 1) ))>. 3 $) 

read(5,400)nchrs, (rlstr(i:i) , i-1,10)  !!  get  string 

format (q, lOal) 

if  (nchrs  .eq.  0)  return  !!  take  default  if  no 
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! !  entry 

call  s2r(rl , rlstr)  !!  convert  str-->int 

if  ((rl  .le.  mx)  .and.  (rl  .ge.  mn) )  then  !!  check  range 

rel-rl 
return 

endif 

lmx-5+max(l , int(loglO(abs (mx)+. 1) ) ) 
lmn-5+max(l , int( loglO(abs(mn)+. 1) ) ) 

write(6 , 300)mn,mx  !!  not  in  range 

00  formatC  The  (min, max) are  ( ' , f<lmn> . 3 , ' , ' , f<lmx>. 3 , 1 ) 
x  -  -  reenter: ' , $) 
go  to  200 
end 


subroutine  s2r(r,str) 

converts  string  --  str  --  to  real  --  r. 

implicit  none 
integer  a.l.i.j 
real*4  r, temp.pt 
char ac ter* (*)  str 

temp— 0 . 
pt-1. 

l-len(str) 

do  100  i— 0,1-1 
a-ichar(str(l-i: 1-i) ) 
if  (a  .eq.  45)  go  to  300 
if  (a  .eq.  48)  go  to  200 
if  (a  .eq.  46)  then 
temp-temp/pt 
pt-1. 
go  to  100 

endif 

if  ((a  .It.  49)  .or.  (a  .gt.  57))  go  to  100 
temp- temp+p  t*float(a-48) 

00  pt-pt*10. 

00  continue 

r-temp 

return 

00  r--temp 


! !  itr-sum{str(i)*10**i) 

! !  a-'-’ 

!  !  a-' O' 

! !  a-' . ' 


! !  a  not  in  1-  -9 


return 

end 


subroutine  clrscrn 
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clears  the  screen 

byte  sbuf (4) 

sbuf (l)-27  !  ESC 
sbuf (2)-91  !  [ 
sbuf (3)-50  !  2 
sbuf(4)-74  !  J 
write(6 , 100) (sbuf (i) ,i-l,4) 
00  format (lh+,4al) 


return 

end 


logical  function  gtyes(pt) 

prompts  user  w/  pt,  returns  .true,  if  'yes' 


implicit  none 
character  ystr 
character*(*)  pt 

gtyes-. false. 
write(6,100)  pt 
)0  format( '  ' ,a,$) 

read(5,200)  ystr 
)0  format (al) 

if  ((ystr  .eq.  'y')  .or.  (ystr 

return 

end 


subroutine  prterr(msg, ios) 

prints  runtime  error  codes 

character*(*)  msg 
integer  ios 

write(6,100)  msg, ios 
)0  format ( '  ' , a, ' : ' , i4) 

return 


! !  write  prompt 
! !  get  response 

.eq.'Y'))  gtyes-. true. 

! !  true  if  first 
! !  character  is  y  or  Y 


! !  write  message  and  code 
!  !  to  terminal 


p 


ft 


i 


% 


« 


s 

•> 
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integer*4 

integer*4 

integer*4 

character*20 

character*80 

integer*4 

logical 


c 


numb  e  r  _o  f  _p  o  int  s 
output_f i le_numb e r 
starting_point 
temp_string 
oname , iname , tfnam 
och(4) , ich 
gtyes 


initialization 


c 


do  i-1 , 4 
och(i)-14+i 
enddo 
ich— 2 

iname— ' data : rxy tes t . r4 ' 
oname- ' data : rxygm . r4 ' 
tfnam-' data: track. fil' 
number_of_points-100 
sl-1 

nbins— 512 
sbin-1 
scl-1 . 
thr-2 . 
pet-. 7 


input 
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100 


* 

Jk  ' 


c 


c 

,r 


call  gtstr(' Enter  the 
call  gtstr(’ Enter  the 
call  gtstr(' Enter  the 
call  gtintC Enter  the 
call  gtint(' Enter  the 
call  gtint(' Enter  the 
call  gtintC Enter  the 
call  gtrelC  Enter  the 
call  gtrelC  Enter  the 
call  gtrel(' Enter  the 
call  otstr( '  ' , 2) 
if  ( .not. gtyes (' Input 
call  otstr ( '  ' , 5) 


output  gram  file  name’ , oname , 80) 
output  track  file  name tfnam, 80) 
input  file  name ' , iname ,80) 
number  of  epochs' ,number_of_points, 700,1) 
starting  epoch  number ’, si , 9999 , 1) 
number  of  fft  bins ', nbins , 99999 , 1) 
starting  bin  number' , sbin, 99999 , 1) 
input  scale  factor' , scl , 99999 .,. 00001) 
pw  threshold' , thr, 99 ., 1. ) 
pw  percent' , pet , 1. ,0 . ) 

s  ok  (y/n) :  '))  go  to  100 


call  opn (och, oname, 'new' ) 
call  opn( ich, iname ,' old' ) 
call  adecini(number_of_points) 


body 


starting_point  -  1 
call  lfcr 
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call  percom  (  0,  number_of_points ,  5) 
do  i  -  1 ,  number_of_points 

call  percom  (  i,  number_of_points ,  5) 
c 
c 

call  adecrd  (  ich, i , ff t_amplitudes , si , scl ,nbins , sbin) 
c 

c  ! !  read  in  an  FFT  line  and  normalize 

c 

starting_point  -  starting_point  +  512 

call  doadec  (  fft_amplitudes)  !  process  FFT  line 

if  (  mod  (  i,  175)  .eq.  1  .and.  i  .ne.  1)  then 
call  abelfin  (  och,gram_number,  gram) 
do  j  -  1,  89600 
gram  (  j)  -  0 
end  do 

gram_number  -  gram_number  +  1 

temp_string  -  'adec*  //  char  (  gram_number  +  48)  //  '  -  grin' 
call  opn(och(gram_number) , temp_string, 'new' ) 

end  if 

call  adecgram  (  mod  (  i  -  1,  175)  +  1,  gram)  !  output  adec  results 
call  writestuff (i , fft_amplitudes , thr ,pct,number_of_points , tfnam) 
end  do 

call  abelfin  (  och,  gram_number,  gram) 
call  els (Ich) 
end 
c 

include  'jsaio.for'  !!  get  jsa's  i/o 

c 
c 
c 
c 
c 

subroutine  abelfin(ch , gn,x) 
c 
c 

c  writes  gram  --  x  to  channel  ch(gn) 

c 

c 

implicit  none 

integer*4  ch(4),gn,i,j 

real  temp(512) ,x(512 , 175) 

c 
c 
c 

do  i-1,512 
temp( i)-0 
enddo 
c 

c  150 
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do  2000  j -1,512 
temp(j)-x(j , i) 
continue 

call  wrec(ch(gn) , temp,4*(i-l)+l,4*i) 
continue 


call  cls(ch(gn)) 
end 

subroutine  assoc  (  association_threshold,  association_gate_min, 

1  association_gate_max,  association_gate_slope ,  fft_amplitude , 

1  fft_associated,  max_ff t_bin_number) 


implicit  none 

real 

assoc iation_gate_max 

real 

association_gate_min  (  1) 

real 

association_gate_slope  (  1) 

real 

association_threshold 

real 

fft_amplitude  (  1) 

logical*l 

fft_associated  (  1) 

integer 

max_f  f  t_b in_numb e  r 

real 

association_gate 

integer 

bin 

logical*l 

lower_edt 

real 

lower_limit 

integer 

lower_limit_bin 

integer 

other_track 

integer 

track 

logical*l 

upper_edt 

real 

upper_limit 

integer 

upper_limit_bin 

include 


'  tracks . inc ' 
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do  bin  -  1,  max_ff t_bin_number 

fft_associaced  (  bin)  -  .false, 
end  do 

do  track  -  1,  number_of_tracks 

if  (  t_adaptive_amplitude  (  track)  .gt.  association_threshold)  then 
association_gate  - 

1  association_gate_min  (  t_type  (  track))  + 

1  association_gate_slope  (  t_type  (  track))  * 

1  (  t_adaptive_amplitude  (  track)  -  association_threshold) 

association_gate  - 

1  min  (  association_gate ,  association_gate_max) 


association_gate  -  association_gate_min  (  t_type  (  track)) 
end  if 


upper_limit  —  min  (  t_frequency  (  track)  +  association_gate , 
float  (  max_ff t_bin_number) ) 

lower_limit  -  max  (  t_frequency  (  track)  -  association_gate ,  1.) 

Check  for  to  see  if  any  higher  frequency  edt  track  will  affect  the 
upper  limit  of  the  association. 

do  other_track  -  track  +  1,  number_of_tracks 
if  (  t_type  (  other_track)  .eq.  edt)  then 
upper_edt  -  .true, 
goto  1000 
end  if 
end  do 

upper_edt  -  . false . 
continue 

if  (  upper_edt)  then 

upper_limit_bin  -  min  (  upper_limit, 

(  t_frequency  (  track)  +  t_frequency  (  other_track) )  *  .5) 
+  0.5 

else 

upper_limit_bin  -  upper_limit  +  .5 
end  if 


Check  to  see  if  any  lower  frequency  edt  track  will  affect 
the  lower  limit  of  the  association 

do  other_track  -  track  -  1,  1,  -1 

if  (  t_type  (  other_track)  .eq.  edt)  then 
lower_edt  -  . true . 
goto  2000 
end  if 
end  do 
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lower  edt  -  .false. 


continue 


if  (  lower_edt)  then 

lower_limit_bin  -  max  (  lower_limit, 

1  (  t_frequency  (  track)  +  t_frequency  (  other  track))  *  .5) 

1  +  .5 

else 

lower_limit_bin  -  lower_limit  +  .5 
end  if 
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Associate  cell(s)  with  track. 

t_associated_amplitude  (  track)  -  fft_amplitude  (  lower_limit_bin) 
t_associated_cell  (  track)  -  lower_limit_bin 
do  bin  -  lower_limit_bin,  upper_limit_bin 

if  (  ff t_amplitude  (  bin)  .gt.  t_associated_amplitude  (  track)) 

1  then 

t_associated_amplitude  (  track)  -  fft_amplitude  (  bin) 
t_associated_cell  (  track)  -  bin 

end  if 

if  (  t_type  (  track)  .eq.  edt)  fft_associated  (  bin)  -  .true. 

!  all  tracks  in  window  of  edt  are  associated 

end  do 

fft_associated  (  t_associated_cell  (  track))  -  .true. 

if  (  t_associated_amplitude  (  track)  .It. 

1  t_adaptive_amplitude  (  track))  then 

t_dcell  (  track)  -  (  t_associated_cell  (  track)  - 
1  t_frequency  (  track))  *  (  t_associated_amplitude  (  track)  / 

1  t_adaptive_amplitude  (  track)) 

else 

t_dcell  (  track)  -  t_associated__cell  (  track)  - 
1  t_frequency  (  track) 


end  if 


end  do 
end 

real  function  besseliO  (  x) 


real  x 

real  arg 

real  b 

integer  i 

real  term 


b  -  1 

arg  -  x  *  x  *  .25 
term  -  1 


do  i  -  1,  25 


P 


term  -  term  *  arg  /  float  (  i)  **  2 

if  (  abs  (  term  /  b)  .It.  le-4)  goto  1000 

b  -  b  +  term 

end  do  , , 
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besseliO  —  b 
end 

subroutine  cmlne  (  window_size,  upper,  lower,  ns,  rx,  x) 

implicit  complex*8  (  a  -  z) 

integer*4  window_size 

integer*4  upper 

integer*4  lower 

integer*4  ns 

real  rx  (  *) 

real  x  (  *) 

integer*4  i 

integer*4  j 

real  list  (  64) 

real  mean 

real  temp 


Build  first  list 

do  i  -  1,  window_size 
list  (  i)  -  rx  (  i) 
end  do 

Sort  list  (  crudely  right  now.) 

do  i  -  1,  window_size 

do  j  -  i  +  1,  window_size 

if  (  list  (  i)  .gt.  list  (  j))  then 
temp  -  list  (  j) 
list  (  j)  -  list  (  i) 
list  (  i)  -  temp 
end  if 
end  do 
end  do 


ft 


v, 

tiS 


v. 

V. 


I  » 

■ 


Calculate  mean 
mean  -  0 . 

do  i  -  lower,  upper 

mean  -  mean  +  list  (  i) 
end  do 

mean  -  mean  /  float  (  upper  -  lower  +  1) 

Calculate  noise  estimate  for  left  end  of  data 

do  i  -  1,  window_size  /  2 
x  (  i)  -  mean 
end  do 

Calculate  noise  estimate  for  middle  of  data 

do  i  -  window_size  /  2  +  1,  ns  -  window_size  /  2 

call  sorto  (  list,  window_size,  rx  (  i  -  window_size  /  2)) 
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call  sorti  (  list,  window_size,  rx  (  i  +  window_size  /  2)) 
call  aver  (  list,  upper,  lower,  mean) 
x  (  i)  -  mean 
end  do 

Calculate  noise  estimate  for  right  end  of  data 

do  i  -  ns  -  window_size  /  2  +  1,  ns 
x  (  i)  -  mean 
end  do 
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include 


tracks . inc ' 


data  free_entry  /  1/ 
data  number  of  tracks  /  0/ 


subroutine  adecfin  (  output_file_number ,  gram) 


implicit 

real 

integer*4 

real 

integer 

integer 

integer 

integer*4 


gram  (  58100) 
output_file_number 

buffer  (  512) 
buffer_pointer 


page  -  0 

buffer_pointer  -  1 

do  i  -  1,  332 

do  j  -  1,  175 

buffer  (  buffer_pointer)  -  gram  ((  j  -  1)  *  332  +  i) 
if  (  buffer_po inter  .eq.  512)  then 

call  ptpage  (  page,  output_f ile_number ,  buffer) 
page  -  page  +  1 
buffer_pointer  -  0 
end  if 

buffer_pointer  -  buffer_pointer  +  1 

end  do 
end  do 

Zero  pad  last  page. 

do  i  -  buffer_pointer ,  512 
buffer  (  i)  -  0. 
end  do 

call  ptpage  (  page,  output_f ile_number ,  buffer) 
call  nclose  (  output_f ile_number ) 


subroutine  adecgrara  (  i,  gram) 


implicit 

real 

integer 

integer 

integer 


gram  (  512,  175) 
i 

bin 

track 
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include  ' tracks. inc' 

do  bin  -  1,  512 

gram  (  bin,  i)  -  0. 
end  do 

do  track  -  1 ,  number_of_tracks 

bin  -  t_frequency  (  track)  +  .5 
if  (  t_detection_flag  (  track) )  then 
gram  (  bin,  i)  -  7. 

else 

gram  (  bin,  i)  -  2. 
end  if 
end  do 


end 

subroutine  adecini(number_of_points) 


implicit  none 


integer*4 

character*80 

logical*l 

integer 

real 

real 

real 

character*80 

character*80 

integer*4 

integer*4 

real 


numb e r _o f _p o i n t s 
coeff icient_file_name 
error 
i 

new_tr ack_p  fa 

old_track_pd 

old_track_pfa 

output_file_narae 

param_file_name 

param_unit 

range  (  2) 

temp_real 


include 

include 


jhj lib . inc ' 
par am. inc ' 
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call  getjfn  (  param_unit) 

„  error  -  .false. 

^.00  param_f ile_narae  -  'adec.prm' 

I  if  (  error)  call  outmes  (  'Couln''t  find  file.  Try  again.') 

error  -  . true . 

■v,  call  askstr  (  'Adec  parameter  file' ,  param_f ile_name) 

open  (  unit  -  param_unit,  err  -  100,  type  -  'old',  readonly, 

1  form  -  'formatted',  file  -  param_file_name) 

P|:  Read  in  adec  paramters 

type  1000 

w  read  (  param_unit,  *)  amplitude_smoothing_constant 

type  1001,  ' amplitude_smoothing_constant' ,  amplitude_smoothing_constant 

read  (  param_unit,  *)  association_gate_min  (  1) 

type  1001,  ' association_gate_min' ,  assoc iation_gate_min  (  1) 

read  (  param_unit,  *)  association_gate_max 

type  1001,  ' association_gate_max' ,  association_gate_max 

do  i  -  1,  3 

.£■  read  (  param_unit,  *)  association_gate_slope  (  i) 

type  1001,  ' association_gate_slope  '  //  char  (i  +  48), 

1  association_gate_slope  (  i) 

'■*  end  do 

read  (  param_unit,  *)  association_threshold 
v  type  1001,  ' association_threshold' ,  association_threshold 

g§  read  (  param_unit,  *)  old_track_pd 

type  1001,  ' old_track_pd' ,  old_track_pd 
, .  read  (  param_unit,  *)  old_track_pfa 

*  type  1001,  ' old_track_pfa' ,  old_track_pfa 

‘‘  detection_threshold  -  log  (  old_track_pd  /  old_track_pfa) 

drop_threshold  -  log  ((1.  -  old_track_pd)  /  (  1.  -  old_track_pfa) ) 
ttj  read  (  param_unit,  *)  temp_real 

display_threshold  -  detection_threshold  +  temp_real 
type  1001,  ' display_threshold' ,  display_threshold 
read  (  param_unit,  *)  dynamic_tracking_threshold 
y!  type  1001,  ' dynamic_tracking_threshold' ,  dynamic_tracking_threshold 

“  read  (  param_unit,  *)  faca 

type  1001,  'faca',  faca 

P  read  (  param_unit,  *)  log_likelihood_max 

type  1001,  ' log_likelihood_max' ,  log_likelihood_max 
read  (  param_unit,  *)  max_fft_bin_nuraber 
type  1002,  'max_fft_bin_number '  ,  max_fft__bin_number 
“v  read  (  param_unit,  *)  merge_gate 

type  1001,  'merge_gate' ,  merge_gate 
read  (  param_unit,  *)  rate_gate 

V  type  1001,  'rate_gate',  rate_gate 

V  read  (  param_unit,  *)  new_track_pfa 

type  1001,  ' new_track_pfa' ,  new_track_pfa 

4  ", 

close  (  unit  -  param_unit) 
type  1003 

$ 

>  do  i  -  1,  3 

new_track_threshold  (  i)  -  sqrt  (  -2.  *  log  (  new_track_pfa) ) 

WV  end  do  . _ . 
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association_gate_min  (  2)  -  association_gate_min  (  1) 
associatlon_gate_min  (  3)  -  association_gate_min  (  1) 

detection_chreshold  -  log  (  old_track_pd  /  old_track_pfa) 
drop_threshold  -  log  ((1.  -  old_track_pd)  /  (  1.  -  old_track_pfa) ) 


a 

a 


000  format  (  ' OAdec  parameters:'//) 

001  format  (  lx,  a,  gl3.6) 

002  format  (  lx,  a,  ,  il3) 

003  format  (  '0') 

end 

subroutine  adecrd  (  ch,  starting_point , 

1  fft_amplitudes ,  si,  scl .nbins , sbin) 


real 

ff t_amplitudes  (  *),scl 

integer*4 

ch , s 1 , i , j , k , nb ins , sb in 

integer*4 

starting_point 

complex 

coefficient_array  (  512) 

real 

nse  (  512) 

real 

temp  (  512) 

real 

temp 2  (  9999) 

include 

' jhj lib . inc' 

r- 


i 

9 


j - ( s 1+s  tar t ing_po int -  2 ) *nb ins+sb in 
k-int(float(j)/128)+l 
1-k+int ( f loat (nb ins ) /128 ) +1 
call  rrec(ch, temp2 ,k, 1) 

ii-j - (k-l)*128+l 

do  i-1,512 

temp(i)-scl*temp2(i+ii) 

enddo 


call  cmlne  (  32,  16,  1,  512,  temp,  nse) 


do  i  -  1,  512 


fft_amplitudes  (  i)  -  scl*temp2  (  i+ii) 


end  do 

end 

program  adecw 

implicit  none 

integer*4 

coeff icient_file_nuraber 

real 

fft_amplitudes  (  512),  scl, 

real 

gram  (  89600) 

integer 

gram  number  /  1/ 

integer 

i.J 

integer*4 

si , sbin, nbins 

/  sqrt(nse(i) )*. 558 


thr,  pet 


160 


a 
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3 

3 

$ 


a 
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subroutine  sorto  (  list,  window_size,  value) 

real  list  (  1) 

integer*4  window_size 

real  value 

integer*4  i 

integer*4  j 

do  i  -  1,  window_size 

if  (  value  ,eq.  list  (  i))  goto  100 
end  do 

00  continue 

Compress  list 

do  j  -  i,  window_size  -  1 

list  (  j)  -  list  (  j  +  1) 
end  do 

end 


vs 
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subroutine  sorti  (  list,  window_size,  value) 


real 

list  (  1) 

integer*^ 

window  size 

real 

value 

integer*4 

i 

integer*4 

j 

real 

tempi 

real 

temp2 

do  i  -  1 , 

window_size  -  1 

if  (  value  .It.  list  (  i))  then 
goto  100 
end  if 
end  do 

i  -  window_size 

.00  continue 

temp2  -  value 

Compress  list 

do  j  -  i,  window_size 
tempi  -  list  (  j) 
list  (  j)  -  temp2 
temp 2  -  tempi 
end  do 

end 


j ii 

M 


a 


P 


a 
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a 
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list  (  1) 

upper 

lower 


real 

integer*4 

integer*4 

real 

mean  -  0. 


do  i  -  lower,  upper 

mean  -  mean  +  list  (  i) 
end  do 

mean  -  mean  /  float  (  upper  -  lower  +1) 
end 

subroutine  compress 


implicit  none 

integer 

integer 

integer 

include 


targe t_track 
track 

trackcounter 
' tracks . inc ' 


track_counter  -  0 
targe t_track  -  1 

do  track  -  1,  free_entry  -  1 

if  (  t_frequency  (  track)  .ge.  0)  then 

track_counter  -  track_counter  +  1 
if  (  track  .ne.  target_track)  then 

t_adaptive_amplitude  (  target_track)  - 
1  t_adaptive_amplitude  (  track) 

t_age  (  target_track)  -  t_age  (  track) 
t_associated_amplitude  (  target_track)  - 
1  t_associated_amplitude  (  track) 

t_associated_cell  (  target_track)  - 
1  t_associated_cell  (  track) 

t_dcell  (  target_track)  -  t_dcell  (  track) 
t_detection_flag  (  target_track)  - 
1  t_detection_flag  (  track) 

t_frequency  (  target_track)  -  t_frequency  (  track) 
t_frequency_rate  (  targe t_track)  - 
1  t_frequency_rate  (  track) 

t_id  (  target_track)  -  t_id  (  track) 

t_log_likelihood  (  target_track)  - 
1  t_log_likelihood  (  track) 

t_type  (  target_track)  -  t_type  (  track) 

end  if 


target_track  -  target_track  +  1 
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end  if 
end  do 

number_of_tracks  -  track_counter 
free_entry  -  number_of_tracks  +  1 


subroutine  deltra  (  track) 


implicit  none 
integer 


track 


include 


' tracks . inc 1 


t_frequency  (  track)  -  -1.0 


subroutine  doade^  (  fft_amplitudes) 
implicit  none 


logical*l 

include 

include 


fft_amplitudes  (  *) 

alpha  (  32)  /  10  *  .18750,  .20099,  .21301,  .22601, 
.23700,  .24899,  .26099,  .28125,  .30600,  .32999, 

.36499,  12  *  .37900/ 

beta  (  32)  /  9  *  .87500,  .01941,  .02246,  .02637, 

.02863,  .03186,  .03640,  .03918,  .04602,  .05527, 

.06522,  .07661,  12  *  .08862/ 

fft_associated  (  512) 

' param . inc ' 

' tracks . inc ' 
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s 

c  Associate  existing  tracks  with  FFT  cells. 


i 


call  assoc  (  association_threshold,  association_gate_min, 

1  association_gate_max,  association_gate_slope ,  fft_amplitudes , 

1  fft_associated,  max_fft_bin_number) 

Update  existing  tracks. 


call  update  (  amplitude_smoothing_constant ,  detection_threshold, 
1  display_threshold,  drop_threshold,  log_likelihood_max) 


c  Create  new  tracks . 

» r. 

.v 

.*»■  call  newtra  (  fft_amplitudes ,  f ft_associated,  max_fft_bin_number , 

1  new_track_threshold) 


5 


% 


Put  the  track  table  back  into  order  of  increasing  frequency, 
call  sort 

Merge  equivalent  tracks . 


-  a. 


call  merge  (  merge_gate,  rate_gate) 

Smooth  and  predict  new  frequency  and  rate. 


call  smooth  (  dynamic__tracking_threshold,  faca,  alpha,  beta, 
1  max_fft_bin_number) 


/- 


Pr 


end 

subroutine  domerg  (  trackl,  track2) 

implicit  none 

integer  trackl 

integer  track2 

integer  deleted_track 

integer  retained_track 


include 


'  tracks . inc ' 
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if  ((  t_type  (  trackl)  .eq.  edt)  .and.  (  t_type  (  track2)  .eq.  edt)) 
1  Chen 

if  (  t_detection_flag  (  trackl)  .eq.  t_detection_flag  (  track2)) 
1  then 

if  (  t_adaptive_ampli tude  (  trackl)  .ge. 

1  t_adaptive_amplitude  (  track2))  then 

retained_track  -  trackl 

else 

retained_track  -  track2 
end  if 

else 

if  (  t_detection_flag  (  trackl))  then 
retained_track  -  trackl 

else 

retained_track  -  track2 
end  if 
end  if 

else  if  (  t_type  (  trackl)  .eq.  edt)  then 
retained_track  -  trackl 
else  if  (  t_type  (  track2)  .eq.  edt)  then 
retained_track  -  track2 

else 

if  (  t_log_likelihood  (  trackl)  .gt.  t_log_likelihood  (  track2)) 
1  then 

retained_track  -  trackl 

else  if  (  t_log_likelihood  (  trackl)  .It. 

1  t_log_likelihood  (  track2))  then 

retained_track  -  track2 

else  if  (  t_adaptive_amplitude  (  trackl)  ,gt. 

1  t_adaptive_amplitude  (  track2))  then 

retained  track  -  trackl 


else 
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1 

s 

i? 

“* '■  m 


end  if 

if  (  retained_track  .eq.  trackl)  then 
deleted_track  -  track2 

else 

deleted_track  -  trackl 
end  if 

if  (  t_age  (  trackl)  .  ge.  t_age  (  track2))  then 
t_id  (  retained_track)  -  t_id  (  trackl) 

else 

t_id  (  retained_track)  -  t_id  (  track2) 
end  if 


c 


.c 

k  a 


c 

i 


N 

.S 

c 

I 


t_age  (  retained_track)  -  max  (  t_age  (  trackl),  t_age  (  track2)) 

call  deltra  (  deleted_track) 

end 

file  i/o  subroutines  -  -  records 
written  as  512  byte  blocks 


subroutine  opn(ch,nam,oon) 

opens  file  w/  name-nam  on  unit  ch 

implicit  none 
integer  ch.ios 
character*^*)  nam 
character*(3)  oon 

! !  open  file 

open ( uni t-ch, name-nam, status-oon, err-100, iostat-ios , 
x  recordtype-' fixed' , recordsize-512/4 , access-' direct ' , 

x  organization-' sequential ' ,blocksize-512*64) 

return 

call  prterr('opn  error’.ios)  !!  print  error 

return 

end 


subroutine  cls(ch) 
closes  file  on  unit  ch 
implicit  none 


‘.'On  ‘ 
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r  • 


C 


c 

100 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


integer  ch.ios 

close (uni t-ch, status-' keep' ,err-100)  !!  close  file 

return 

call  prterr('cls  error', ios)  !!  print  error 

return 

end 


subroutine  rrecxy(chx,xbuf , chy ,ybuf , sr , fr ,xey) 


reads  starting  sr  finnishing  fr:  chx-->xbuf,  chy-->ybuf 
if  xey  sets  ybuf-xbuf  --no  chy  exsists 


implicit  none 

integer  chx, chy ,nrcds , sr , fr , i 

byte  xbuf (abs( (fr-sr+l)*512) ) ,ybuf (abs((fr-sr+l)*512) ) 
logical  xey 


b' 


$ 


5 


l 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


endif 

return 

end 


subroutine  rrec(ch,buf , sr , fr) 

reads  fr-sr  512  byte  records  from  ch  to  buf,  starting  at 
record  sr  finishing  at  fr 

implicit  none 

integer  sr , fr , ch, ios , i , s , j 

byte  buf (abs((fr-sr+l)*512) ) , temp(512) 


if  (sr  .gt.  fr)  then 
ios-999 
goto  200 

endif 


! !  starting  record  after 
! !  ending  record 
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u  . 

call  rrec(chx, xbuf , sr.fr) 

!!  read  x's  records 

ft 

y  c 

if  (xey)  then 

U  xey-->set  y-x 

d 

nrcds-fr-sr+1 

■J; 

do  200  i-l,nrcds*512 

r,* 

.',r 

:::  200 

ybuf ( i ) -xbuf ( i ) 

r 

# 

else 

s . 

& 

call  rrec(chy,ybuf , sr.fr) 

!!  else  read  y's  records 

& 


■■V 

<V 


8 


& 


■c. 


i* 


* 
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!!  read  records  sr...fr 


I 


read(unit-ch , rec-i , err-200 , iostat-ios ) temp 


s-(i-sr)*512 
f-.’  do  j -1,512 

buf  ( j+s)-temp(  j  ) 
enddo 


put  read  buffer  in 
return  variable,  buf 
512  bytes/record 


continue 

return 


* 

t«  v 

'•'c 

c 


£ 


I 


£ 


p 


if0 

i 

200 
K  j 


C 


call  prterr(’rrec  err’, ios)  !!  print  errors 

return 

end 


subroutine  wrec(ch,buf ,sr , fr) 

writes  buf  to  ch  in  512  byte  records  -- 
starting  at  record  sr,  finnishing  at  fr 

implicit  none 

integer  s ,ch,sr , fr , ios , i, j 

byte  buf (abs( (fr-sr+l)*512) ) , temp(512) 


if  (sr  .gt.  fr)  then 
ios-999 
go  to  200 

endif 

do  100  i-sr,fr 

s-(i-sr)*512 
do  j-1,512 
temp( j )-buf (j+s) 
enddo 


!  !  starting  record  after 
!  !  last  record 


!!  write  records  sr...fr 


!!  fill  write  buffer  -- 
!  !  512  bytes/record 


write (uni t-ch , rec-i , err-200 . iostat-ios ) temp 
return 


call  prterr(’wrec  error’, ios)  !!  print  errors 

return 

end 


c  terminal  i/o  routines 

« 


169 


7*7. 


DRA2 : [ DANIEL] HP . TMP ; 1 


24-FEB-1986  12:28 


subroutine  otstr(str ,nlf) 

prints  str,  on  the  current  line,  followed  by  nlf  cr/linefeeds 


implicit  none 
integer  nlf,l 
character*(*)  str 

write(6 , 100)str  !!  prints  string  to  terminal 

format(lh+,a, $) 

if  (nlf  .It.  1)  return  !!  prints  If's 

do  200  1-1, nlf 

write(6 , 300) 

format( 1  ' ) 

continue 

return 

end 


subroutine  gtstr(pt,var ,n) 

gets  a  string  from  the  terminal  --  prompts  w/  pt,  var-string 

var  has  length  n<101 

implicit  none 
character*(*)  pt.var 
integer  n,i,nchrs 
character*100  invar 

i-min(len(var) , index(var , '  ')) 

write(6,100)  pt.var  ! !  write  prompt,  default 

formate  '.a,'  ( ' ,a<i>, ' ) : ' ,$) 

read(5 , 200) nchrs ,( invar (i : i) , i-1 ,nchrs)  !!  read  srting 
format (q , lOOal) 

if  (nchrs  ,eq.  0)  return  !!  return  default 

do  300  i— 1 , n 
var ( I : i ) - '  ' 

if  (i  . le .  nchrs)  var (i : i)-invar ( i : i)  !!  fill  return  variable 

return 

end 


subroutine  gtint(pt, igr ,mx,mn) 

gets  integer  from  terminal  --  prompts  user  w/  pt,  accepts 

igr  in  range  (mn,mx) 


* 

.  i  . 


a 


U 


V'] 

V'. 


,r> 

s'. 

v> 


C 
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implicit  none 
character*10  itstr 

integer  igr , mx , mn , it , nchrs , i , lmn , lmx 
character*(*)pt 

i-2+max(l , int(loglO(abs (float ( igr) )+ . 1) ) ) 
write(6 , 100)pt , igr  ! 

00  formate  \a,'  ( ’  ,  i<i> 

00  read(5 ,400)nchrs , (itstr(i : i) , i-1 , 10)  ! 

00  format(q, lOal) 

if  (nchrs  .eq.  0)  return  ! 

! 

call  s2i(it, itstr)  ! 

if  ((it  .le.  mx)  .and.  (it  . ge .  mn) )  then 
igr-it 
return 

endif 

lmn-2+raax(l , int (logl0(abs( float j (mn) )+. 1) ) ) 
lmx-2+max(l , int (logl0(abs (float j (mx) )+. 1) ) ) 

write(6 , 300)mn,mx  !!  not  in  range 

00  formatC  The  (min, max) are  (' , i<lmn>, ' , ' , i<lmx>, ' )  --  reenter :',$) 
go  to  200 
end 


subroutine  s2i(itr,str) 

converts  string  -  -  str  -  -  to  integer  -  -  itr . 

I  implicit  none 

integer  a,l 

integer  itr , i ,j , temp ,pt 
character*(*)  str 

temp-0 

pt-1 

i  l-len(str) 

i  do  100  i-0,1-1  !!  itr-sum( str(i)*10**i } 

|  a-ichar(str(l- i : 1- i) ) 

if  (a  .eq.  48)  go  to  200  !!  a-'0’ 

if  (a  .eq.  45)  go  to  300  !!  a-'-' 

if  ((a  .It.  49)  .or.  (a  .gt.  57))  go  to  100  !!  a  not  in  1--9 

temp- temp+p t* ( a - 4 8 ) 

00  pt-pt*10 

00  continue 

I  itr-temp 

1  return 


write  prompt,  default 


get  string 

take  default  if  no 
entry- 

convert  str-->int 

! !  check  range 


00 


itr-- temp 
return 
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subroutine  gtrel(pt , rel , mx,mn) 

gets  real*4  from  terminal  --  prompts  user  w/  pt,  accepts 

rel  in  range  (mn,mx) 


implicit  none 
character*10  rlstr 
real*4  rel,mx,mn,rl 
integer  nchrs , i , lmx , lmn 
character* (*)pt 

write(6 , 100)pt , rel  !!  write  prompt,  default 

formate  ',a,’  ( ' , f<5+max(l, int(loglO(abs(rel)+ . 1) ) )>. 3 , ' ) : ' ,$) 

read (5 ,400) nchrs , (rlstr(i : i) , i— 1 , 10)  !!  get  string 

format (q, 10a!) 


if  (nchrs  .eq.  0)  return 
call  s2r(rl , rlstr) 


!  take  default  if  no 
!  entry 

!  convert  str-->int 


if  ((rl  . le .  mx)  .and.  (rl  .ge.  mn))  then  !!  check  range 

rel-rl 
return 

endif 

lmx-5+max(l, int(loglO(abs(mx)+. 1))) 
lmn— 5+raax(l , int (loglO (abs (mn)+ . 1) ) ) 

wri te ( 6 , 300 ) mn , mx  ii  not  in  range 

formate  The  (min.max)are  (' , f<lmn>. 3 , ' , ' , f<lmx>. 3 , ’ ) 

-  -  reenter: ' , $) 
go  to  200 
end 


subroutine  s2r(r,str) 

converts  string  --  st: 

implicit  none 
integer  a,l,i,j 
real*4  r,temp,pt 
character*(*)  str 

temp-0 . 

pt-1. 

l-len(str) 


•  -  to  real  -- 


do  100  i-0,1-1 
a-ichar ( s  tr ( 1 - i : 1- i) ) 
if  (a  .eq.  45)  go  to  300 
if  (a  .eq.  48)  go  to  200 


!!  itr-sum{ str ( i)*10**i ) 


! !  a- ' 0 ' 


‘-A.  /  VaV.V.V 


,  *  .  *’•  " 
-  V.-. 


v  .• -- •  -  ■  -  v.v  •  j' •  -  v 


J-V-Vu' 


y 
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‘200 

100 


I 


poo 

c 


1  V 

& 

c 


V 


V 
•  • 
Vw 


LOO 


'S 


3 

\V 


i 

^00 


#00 

>> 


*1 

i 


if  (a  .eq.  46)  then 
temp-temp/pt 
pt-1. 
go  to  100 

endif 

if  ((a  .It.  49)  .or.  (a  ,gt.  57))  go  to  100  !!  a  not  in  1--9 

temp-temp+pt*float(a-48) 

pt-pt*10 . 

continue 

r-temp 


return 

r--temp 


return 

end 


subroutine  clrscrn 
clears  the  screen 


byte  sbuf(4) 


sbuf(l)-27 
sbuf (2)-91 
sbuf (3)-50 
sbuf (4) -74 


ESC 

[ 

2 

J 


write(6 , 100) (sbuf (i) , i-1 ,4) 
format ( lh+ , 4al ) 


return 

end 


logical  function  gtyes(pt) 

prompts  user  w/  pt,  returns  .true,  if  ’yes’ 


implicit  none 
character  ystr 
character^*)  pt 


! !  write  prompt 
! !  get  response 


gtyes- . false . 
write(6,100)  pt 
formate  *  ,a,$) 
read(5,200)  ystr 
format(al) 

if  ((ystr  .eq.  *y')  .or.  (ystr  .eq.'Y'))  j;tyes- .  true . 

return  t .  true  if  first 

end  !  !  character  is  y  or  Y 
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c 

c 

c 

c 

c 

c 

c 

c 


c 

100 

c 


subroutine  prterr(msg, ios) 

prints  runtime  error  codes 

character*(*)  msg 
integer  ios 

write(6,100)  msg, ios 
format ( '  ' ,a, ' : ' , i4) 

return 

end 


!  !  write  message  and  code 
!  !  to  terminal 


c 


real  function  log_likelihood_ratio  (  type,  amplitude) 


real  amplitude 

integer  type 


real 


besseliO 


Adec  log  likelihood  ratio.  Pfa  -  le-4,  Pd  -  .5 
log_likelihood_ratio  -  log  (  besseliO  (  amplitude))  -  .5 
end 

subroutine  maktra  (  bin,  amplitude,  type) 


implicit  none 


real 

integer 

real 

integer 


amplitude 

bin 

log_likelihood_ratio 

type 


integer  id  /  0/ 


include 


’  tracks . inc ' 


t_frequency  (  free_entry)  -  bin 
t_frequency_rate  (  free_entry)  -  0 . 
t_log_likelihood  (  free_entry)  - 
1  log_likelihood_ratio  (  type,  amplitude) 
t_adaptive_amplitude  (  free_entry)  -  amplitude 
t_age  (  free_entry)  -  0 
t_dcell  (  free_entry)  -  0 
t_type  (  free_entry)  -  type 
t_detection_flag  (  free_entry)  -  0 
t_id  (  free_entry)  -  id 

id  -  id  +  1 

free_entry  -  free_entry  +  1 

end  174 
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subroutine  merge  (  merge_gate,  rate_gate) 
implicit  none 


integer 

integer 

include 


merge_gate 

rate_gate 

other_track 

track 

' tracks . inc 1 


m 


Sn9 


\ « \  -  *r**  *  *  •'»  ^ 

>  V  V  “>  • 
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do  track  -  2 ,  number_of_tracks 

do  other_track  -  track  -  1,  1,  -1 

if  (  t_frequency  (  track)  .It.  0)  goto  1000 

if  (  t_frequency  (  other_track)  .It.  0)  then  !  deleted  track 

else  if  (((  t_type  (  track)  .eq.  weak)  .and. 

1  (  t_type  (  other_track)  .eq.  strong))  .or. 

1  ((  t_type  (  track)  .eq.  strong)  .and. 

1  (  t_type  (  other_track)  .eq.  weak)))  then 

else  if  (  t_frequency  (  track)  -  t_frequency  (  other_track) 

1  .gt.  merge_gate)  then 

goto  1000 

else  if  ((  t_type  (  track)  .eq.  edt)  .and. 

1  (  t_type  (  other_track)  .eq.  edt))  then 

if  ((  t_age  (  track)  .eq.  1)  .or. 

1  (  t_age  (  other_track)  .eq.  1))  then 

call  domerg  (  track,  other_track) 

else  if  (  abs  (  t_frequency_rate  (  track)  - 
1  t_frequency_rate  (  other_track) )  .le.  rate_gate)  then 

call  domerg  (  track,  other_track) 

end  if 

else 

call  domerg  (  track,  other_track) 
end  if 
end  do 

000  continue 

end  do 

call  compress 
end 

subroutine  newtra  (  ff t_amplitude ,  ff t_associated,  max_ff t_bin_number , 
1  threshold) 

implicit  none 

real  fft_amplitude  (  1) 

logical*l  fft_associated  (  1) 

integer  max_ff t_bin_number 

real  threshold  (  1) 


24-FEB- 1986  12:28 

bin 

looking_for_downslope 
1  tracks . inc ' 
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integer 

logical*! 


include 
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looking_for_downslope  -  .false. 

do  bin  -  1,  max_ff t_bin_number  -  1 

if  (  ff t_amplitude  (  bin)  -  fft_amplitude  (  bin  +1)  .gt.  0) 

1  then 

if  (  looking_for_downslope)  then 

if  (  ff t_amplitude  (  bin)  .It.  threshold  (  weak))  then 
else  if  (  fft_associated  (  bin))  then 

else  if  (  fft_amplitude  (  bin)  .gt.  threshold  (  strong)) 
1  then 

call  maktra  (  bin,  fft_amplitude  (  bin) ,  strong) 

else 

call  maktra  (  bin,  fft_amplitude  (  bin),  weak) 
end  if 
end  if 


else 


looking_for_downslope  -  .true, 
end  if 


end  do 


call  compress 


end 

subroutine  smooth  (  dynamic_tracking_threshold,  faca,  alpha,  beta, 
1  max_fft_bin_number) 

implicit  none 


real 

real 

real 

real 

integer 


alpha  (  0:*) 
beta  (  0:*) 

dynamic_tracking_threshold 

faca 

max  fft  bin  number 


integer  i 

integer  track 


include 


tracks . inc ' 


nwvwn  uRUfvrva 


? 

\ 11 
5 

■.* 

|K 
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r 

i 


H 


: 


I 


K 


■* 


5 


s  i 


do  track  -  1 ,  number_of_tracks 

i  -  faca  *  t_adaptive_amplitude  (  track) 

if  (  t  .It.  0)  then 
i  -  0 

else  if  (  i  .gt.  31)  then 
i  -  31 


■ 

end 

if 

if 

(  t_dcell  (  track)  .ne.  0)  then 

t_frequency  (  track)  -  t_frequency  (  track)  + 

,  * 
v  * 

i 

alpha  (  i)  *  t_dcell  (  track) 

if  (  t_adaptive_amplitude  (  track)  .ge. 

m 

i  * 

i 

dynamic_tracking_threshold)  then 

t_frequency_rate  (  track)  -  t_frequency_rate  (  track) 

i 

beta  (  i)  *  t_dcell  (  track) 

end  if 

end 

if 

t_frequency  (track)  -  t_frequency  (  track)  + 

i 

t_frequency_rate  (  track) 

if 

((  t_frequency  (  track)  .gt.  max_fft_bin_number)  .or. 

i 

i 

(t_frequency  (  track)  .It.  1.))  call  deltra  (  track) 

end 

do 

call  compress 
end 

subroutine  sort 
implicit  none 


integer 

integer 

include 


trackl 

track2 

' tracks . inc 1 


do  trackl  -  1,  number_of_tracks  -  1 

do  track2  -  trackl  +  1,  number_of_tracks 

if  (  t_frequency  (  trackl)  .gt.  t_frequency  (  track2))  then 

call  swap_real  (  t_adaptive_amplitude  (  trackl) , 

1  t_adaptive_amplitude  (  track2)) 

call  swap_integer  (  t_age  (  trackl),  t_age  (  track2)) 

call  swap_real  (  t_associated_amplitude  (  trackl) , 

1  t_associated_amplitude  (  track2)) 

call  swap_real  (  t_associated_cell  (  trackl) , 

1  t_associated_cell  (  track2)) 
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9 


call  swap_real  (  t_dcell  (  trackl) ,  t_dcell  (  track2)) 


call  swap_logicall  (  t_detection_flag  (  trackl) , 
t_detection_flag  (  track2)) 


call  swap_real  (  t_frequency  (  trackl) , 
t_frequency  (  track2)) 


call  swap_real  (  t_frequency_rate  (  trackl), 
t_frequency_rate  (  track2)) 


call  swap_integer  (  t_id  (  trackl),  t_id  (  track2)) 


call  swap_real  (  t_log_likelihood  (  trackl) , 
t_log_likellhood  (  track2)) 


call  swap_integer  (  t_type  (  trackl),  t_type  (  track2)) 


end  if 


end  do 


end  do 


m 


vJ,-lQ,%vvTlv^v;t.vA>v>- 
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subroutine  swap_real  (  x,  y) 


■ 


real  x 
real  y 


real  z 
z  -  x 

x  -  y 

y  -  z 


end 

-.t 

iV 


1 


subroutine  swap_integer  (  x,  y) 

integer  x 
integer  y 

integer  z 


z  -  x 

x  -  y 

y  -  z 
end 


subroutine  swap_logicall  (  x,  y) 


I 


* 


% 


logical*l  x 

logical*l  y 

logical*l  z 

z  -  x 

x  -  y 

y  -  z 
end 

subroutine  update  (  amplitude_smoothing_constant,  detection_threshold, 
1  display_threshold,  drop_threshold,  log_likelihood_max) 

implicit  none 


$ 

2 

I 

i 


real 

real 

real 

real 

real 

real 

real 

real 

integer 


amplitude_smoothing_constant 

detection_threshold 

display_threshold 

drop_threshold 

log_likelihood_max 

log_likelihood_ratio 

til 

temp_real 

track 
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do  track  -  1,  number_of_tracks 

t_log_likelihood  (  track)  -  t_log_likelihood  (  track)  + 

1  log_likelihood_ratio  (  t_type  (  track) , 

1  t_associated_amplitude  (  track) ) 

t_log_likelihood  (  track)  - 

1  min  (  t_log_likelihood  (  track)  ,  log_likelihood_max) 


til  -  t_log_likelihood  (  track) 
if  (  til  .le.  drop_threshold)  then 
call  deltra  (  track) 

else  if  (  til  .It.  display_threshold)  then 

t_detection_flag  (  track)  -  .false, 
if  (  t_type  (  track)  .eq.  edt) 

t_age  (  track)  -  t_age  (  track)  +  1 


else  if  (  til  .It.  detection_threshold)  then 
if  (  t_detection_flag  (  track) )  then 

t_age  (  track)  -  t_age  (  track)  +  1 


if  (  t_type  (  track)  .eq.  edt) 

t_age  (  track)  -  t_age  (  track)  +  1 


end  if 


if  (  t_type  (  track)  .eq.  edt)  then 

t_detection_flag  (  track)  -  .true. 
t_age  (  track)  -  t_age  (  track)  +  1 


t_type  (  track)  -  edt 
t_age  (  track)  -  1 
t_detection_flag  (  track)  -  .true. 

end  if 

end  if 

temp_real  -  t_associated_amplitude  (  track)  - 
1  t_adaptive_amplitude  (  track) 
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if  (number_of_tracks . ge .  1)  write(lun, 891)  k 
jjS91  format( ' ( ' , i4) 


main  loop 


do  100  i  -  1,  number  of  tracks 


find  fl , fh, diffusion 


^  fh  first 

j-t_associated_cell(i) 

nmr-0 

noh-1 

nom-0 

776  j-j+1 

if  (fft_amplitudes( j ) . ge . thr*level)  then 
noh-noh+1 
nmr-0 

else 

nom-nom+1 
£  nrar-nmr+1 

if  ( (nmr . ge . 2) . and. (pet . ge . ( float (noh) /float (noh+nom) ) ) ) 

■v  x  go  to  777 

jj  endif 

go  to  776 

777  continue 

m  j-j-i 

{P  if  ( (fft_amplitudes(j ) . It . thr*level) .and. (j .ne . t_associated_cell(i) ) ) 

x  go  to  777 

• .  fh-j 


do  fl 


j-t_associated_cell(i) 

>j,  nmr-0 

V  noh-1 

^  nom— 0 

778  j-j-1 

if  (fft_amplitudes(j  ) .  ge .  thr*level)  then 

V  noh-noh+1 
nmr-0 

else 

^  nom-nom+1 

nmr-nmr+1 

M  if  ((nmr.ge.2) .and. (pct.ge. (float(noh)/float(noh+nom) ) ) ) 

wv  x  go  to  779 

endif 
go  to  778 
0C79  continue 

g  185 


ft 


o  o 


DRA2: [ DANIEL] HP. TMP;1 


24-FEB-1986  12:28 


j-j+1 

if  ( (ff t_amplitudes (j ) . It . thr*level) . and. ( j . ne . t_associated_cell (i) ) ) 
x  go  to  779 

fl-j 
c 
c 

c  find  diffusion 

c 

c 

dif-0 

do  J-fl,fh 

dif-fft_amplitudes(j )+dif 
enddo 

dif-dif/( float (fh-f 1+1)) 


C 

c 

c 

c 


c 

c 

c 

c 

c 

c 

C1001 

c 

c 

c 

c 

c 


889 

c 

c 

c 

c 

c 

100 

c 

890 


type  1001,  i, 

1  t_detection_flag  (  i) , 

1  t_frequency  (  i) , 

1  t_frequency_rate(i) , 

1  t_adaptive_amplitude  (i), 

1  t_associated_cell  (  i) , 

1  t_associated_amplitude  (  i) , 

1  fl. 

1  fh, 

1  dif , 

1  level 

format  (  i4 , 15 , f 8 . 2 , f 8 . 2 , f7 . 2 , f8 . 0 , f8 . 2 , i7 , i7 , f7 . 2 , f 7 . 2) 


write  stuff  to  output  file 


write ( lun, 889)  t_detection_flag  (  i) , 

1  t_frequency  (  i) , 

1  t_frequency_rate( i) , 

1  t_adaptive_amplitude  (i), 

1  t_associated_cell  (  i) , 

1  t_associated_amplitude  (  i) , 

1  fl, 

1  fh, 

1  dif, 

1  leve 1 

format  (  '  ( ' , 15 , f8 . 2 , f8 . 2 , f7 . 2 , f 8 . 0 , f8 . 2 , i7 , i7 , f7 . 2 ,  f7 . 2  ,  ' ) ' ) 


continue 

if  (number_of_tracks . ge . 1)  write (lun, 890) 
format ( ' ) ' ) 


* 
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close  output  file 
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if  (j.eq.nop)  close (unit-lun, status-' keep 1 ) 


m 
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JHJLIB . INC 

- 

integer*4 

parameter 

ap_max_addr e  s  s 
(  ap_max_address  -  16383  ) 

i 

integer*4 

parameter 

array_processor 
(  array_processor  -  2  ) 

> 

V 

1 

integer*4 

parameter 

askflt  1 
(  askf lt_l  -  1  ) 

i 

integer*4 

parameter 

askflt_luv 
(  askflt_luv  -  7  ) 

! 

1 

integer*4 

parameter 

askflt_lv 
(  askflt_lv  -  5  ) 

! 

i 

■ 

integer*4 

parameter 

askflt_lu 
(  askflt_lu  -  3  ) 

i 

1 

integer*4 

parameter 

askf 1 t_node  f aul t 
(  askflt_node fault  -  0  ) 

integer*4 

parameter 

askflt_u 
(  askflt_u  »  2  ) 

» 

j 

integer*4 

parameter 

askflt_uv 
(  askflt_uv  -  6  ) 

1 

integer*4 

parameter 

askflt_v 
(  askflt_v  -  4  ) 

i 

i 

lnteger*4 

parameter 

askint_l 
(  askint_l  -  1  ) 

integer*4 

parameter 

askint_luv 
(  askint_luv  -  7  ) 

integer*4 

parameter 

askint_lv 
(  askint_lv  -  5  ) 

integer*4 

parameter 

askint_lu 
(  askint_lu  -  3  ) 

j 

integer*4 

parameter 

askint_node fault 
(  askint_nodefault  -  0  ) 

integer*4 

parameter 

askint_u 
(  askint_u  -  2  ) 

lnteger*4 

parameter 

askint_uv 
(  askint_uv  -  6  ) 

integer*4 

parameter 

askint_v 
(  askint_v  -  4  ) 

integer*4 

parameter 

comp lex_floa ting 
(  complex_floating  -  3  ) 
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V 

integer*^ 

complex_integer 

i 

parameter 

(  complex_integer  -  4  ) 

integer*4 

memory 

!v 

parameter 

(  memory  -  1  ) 

integer*4 

read_access 

■ 

parameter 

(  read_access  -  1  ) 

integer*4 

real_floating 

parameter 

(  real_floating  -  1  ) 

■s 

integer*4 

real_integer 

parameter 

(  real_integer  -  2  ) 

integer*4 

write_access 

parameter 

(  write_access  -  2  ) 

V, 

real 

amplitude  smoothing  constant 

real 

association_gate_min  (  3) 

real 

association_gate_max 

w  * 

real 

association  gate_slope  (  3) 

V 

real 

associationthreshold 

V 

real 

detectionthreshold 

real 

display_threshold 

«,* 

real 

dropthreshold 

n 

real 

dynamic_tracking_threshold 

real 

faca 

» 1 

real 

log  likelihood  max 

.  •* 

integer 

max_f  f  t_b innumbe  r 

** 

real 

merge_gate 

real 

new_track_threshold  (  3) 

l 

real 

rate_gate 

common  /  param  / 

K 

1  amplitude_smoothing_constant, 

v 

1  association_gate_min. 

I  P* 

1  association_gate_max, 

1  assoc iation_gate_s lope , 

jti 

1  association_threshold, 

o' 

1  detection  threshold, 

1  display 

threshold, 

1  drop  threshold, 

1  dynamic 

_tracking_threshold, 

1  faca, 

,  n- 

1  log  likelihood  max, 

.v# 

1  max_fft_ 

bin  number, 

.1* 

1  merge_gate, 

1  new_track_threshold, 

V*W 

1  • 

1  rate_gate 

Jv 

c 

tracks . inc 

$ 

integer 

edt 

-V 

parameter 

(  edt  -  3) 

integer 

max_tracks 

V, 

parameter 

(  max_tracks  -  400) 

■ 
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integer 

strong 

parameter 

(  strong  -  2) 

integer 

weak 

parameter 

(  weak  -  1) 

integer 

free_entry 

integer 

number_of_ tracks 

real 

t_adaptive_amplitude  ( 

max_tracks) 

integer 

t_age  (  max_tracks) 

real 

t  associated  amplitude 

(  max_tracks) 

real 

t  associated  cell  (  max  tracks) 

real 

t_dcell  (  max_tracks) 

logical*l 

t_detection_flag  (  max_ 

tracks) 

real 

t_frequency  (  max_tracks) 

real 

t_frequency_rate  (  max_ 

tracks) 

integer 

t_id  (  max_tracks) 

real 

t_log_likelihood  (  max_ 

_tracks) 

integer 

t_type  (  max_tracks) 

common  / 

tracks/ 

I  free_entry, 

1  number_of_tracks , 

1  tadap  t ive_amp 1 i tude , 

1  t_age , 

1  t_associated_amplitude , 

1  t_associated_cell , 

1  t_dcell , 

1  t_detection_flag, 

1  t_frequency  , 

1  t_f requency_rate , 

1  t_id, 

1  t_log_likelihood, 

1  t_type 


24-FEB- 1986  12:32 


,v 


DRA2 : [ DANIEL] HP . TMP ; 3 


>3 

'■'B.O 


3.0 

|:o° 

0.5 

^2.25 

L>3.5 

le-4 

1:°5 

4.0 


,-ylO .  017 


$12 

D  .  0 


0.5 

.03 


v.’ 


.> 


fB 


§ 


■> 


*-3 

3 


x  T 


(weak) 

(strong) 

(edt) 


amp 1 i tude_smoo  thing_cons  tant 
association_gate_min 
assoc iation_gate_max 
assoc iation_gate_s lope 
association_gate_slope 
association_gate_slope 
assoc iation_thresho Id 
old_track_pd 
old_track_pfa 

delta_display_threshold_from_detection_threshold 

dynamic_tracking_threshold 

faca 

log_l ike 1 ihood_max 
max_f  f  t_b in_number 
merge_gate 
rate_gate 
new_track_pfa 
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